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Preface 

The Center for Modeling of Turbulence and Transition (CMOTT), a co- 
operative turbulence research team, was formally established in May of 1990 
as a focal group within the Institute of Computational Mechanics for Propul- 
sion (ICOMP). The location of CMOTT is shown in the organization struc- 
ture chart in Appendix A. The objectives of the CMOTT are to develop, vali- 
date and implement the models for turbulence and boundary-layer transition 
in practical engineering flows. The flows of interest are three-dimensional, 
incompressible and compressible flows with chemistry. The schemes being 
studied include the two-equation (e.g. k-e) and algebraic Reynolus-stress 
models, the full Reynolds-stress (or second moment closure) models, the 
probability density function (pdf) models, the Renormalization Group The- 
ory (RNG) and Direct Interaction Approximation (DIA), the Large Eddj' 
Simulation (LES) and Direct Numerical Simulation (DNS). 

Currently, CMOTT has eight formal members working on various as- 
pects of turbulence and transition modeling in collaboration with NASA- 
Lewis scientists and Case Western Reserve University (CWRU) faculty mem- 
bers. The CMOTT members have been actively involved in international and 
national turbulence research activities through meetings, seminars, work- 
shops and exchange-visitors. Since June of 1990, a CMOTT seminar series 
has been conducted with speakers invited from within and outside of the 
NASA Lewis Research Center, including foreign speakers. In 1991, a new 
series of biweekly CMOTT technical meetings was initiated for informal dis- 
cussions regarding special issues in turbulence and transition modeling. The 
CMOTT research activity is advised by a group consisting of Professor J.L. 
Lumley (Cornell University), Dr. M. Goldstein (NASA/LeRC) and Professor 
E. Reshotko (CWRU). 

This research brief contains the progress reports of the CMOTT Re- 
search staff from May 1990 to May 1991. It is intended to be an informational 
report of the CMOTT activities as well as an annual report to ICOMP and 
NASA. The current CMOTT roster and its organization are listed in the 
Appendix A. Listed in Appendix B.l are the visiting members and their 
seminar abstracts. Appendix B.2 gives the scientific and technical issues dis- 
cussed in biweekly CMOTT meetings. Journal and conference publications 
by CMOTT members are grouped in Appendix C. 
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Starting in 1991, NASA Technical Memoranda authored by members of 
the CMOTT staff will be given a specific number to identify them as CMOTT 
reports. These manuscripts will be made available for early dissemination of 
completed research results by the CMOTT staff. 

Finally, we express our thanks to one of the CMOTT members, Dr. 
William W. Liou, who carefully assembled the material and provided edito- 
rial assistance. 


Louis Povinelli 
Meng-Sing Liou 
Tsan-Hsing Shih 
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The Study of PDF Turbulence Models in Combustion 

Andrew T. Hsu 


1. Motivation and Objectives 

1.1 Motivation 

The accurate prediction of turbulent combustion is still beyond reach 
for today’s computation techniques. It is the consensus of the combustion 
profession that the predictions of chemically reacting flow were poor if con- 
ventional turbulence models were used. The main difficulty lies in the fact 
that the reaction rate is highly non-linear, and the use of averaged temper- 
ature, pressure and density produces excessively large errors. The proba- 
bility, density function (pdf) method is the only alternative at present time 
that uses local instant values of the temperature, density, etc., in predicting 
chemical reaction rate, and thus is the only viable approach for turbulent 
combustion calculations. 

1.2 Objectives 

The present work aims at the development and implementation of the 
the pdf turbulence models in solving realistic combustion problems. The fact 
that the pdf equation has a very large dimensionality renders finite difference 
schemes extremely demanding on computer memories and thus impractical, 
if not entirely impossible. A logical alternative is the Monte Carlo scheme, 
which has been used extensively in statistical physics. However, the use of 
Monte Carlo scheme to solve both the fiowfield and the chemical reaction is 
very time consuming. Further more, since CFD has reached a certain degree 
of maturity as well as popularity, it seems less beneficial to abandon CFD 
completely and opt for Monte Carlo schemes. Therefore, we propose the 
use of a combined CFD and Monte Carlo scheme in the present study. The 
scheme would use the conventional flow solvers when calculating the fiowfield 
properties such as velocity, pressure, etc., while the chemical reaction part 
would be solved using Monte Carlo solvers. 
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2. Works Accomplished 

2.1 Code development. 

A parabolic code with k — e turbulence models have been developed 
in the past months. Three different k — e models have been tested with 
satisfactory numerical results. 

A grid dependent Monte Carlo scheme is being explored. This scheme 
discretize the pdf equation on a given grid and write, for parabolic flows: 

Px+dx,j = a jPx,j+i + fijPx,j + IjPxj-l (1) 

and we require 

a j + 0j + 7j = 1 (2) 

Using a very simple test case of a convection/diffusion process with two 
scalars, it was found that the previous scheme does not conserve mass frac- 
tions due to re-contamination. It is found that in order to conserve the 
mass fractions absolutely, one needs to add further restriction to the scheme, 
namely 

aj + 7i = oij-! + 7i+i (3) 

A new algorithm was devised and is currently being tested. Again using 
the simple test case of two scalars with assumed constant coefficients in the 
pdf equation, the new algorithm is shown to conserve the mass fractions 
perfectly in cases of uniform flows or pure diffusion problems. Deficiencies 
such as directional bias and re- contamination that were found in the previous 
algorithm are completely eliminated. 

2.2 Applications. 

The code developed has been validated by solving a heated turbulent 
jet. The temperature is treated as a conserved passive scalar and solved 
using the pdf Monte Carlo simulation while the flow field is obtained using 
a conventional CFD solver. The mean temperature profile and RMS of the 
temperature fluctuation were compared with experimental data. 

As a first application to combustion problems, the non-premixed flame 
of hydrogen and fluorine is being studied. A comparison between primary 
results from the present study and experimental data show that the present 
scheme predicts the mean flame temperature accurately. 
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3. Future Plans 

1. Further investigate the case of hydrogen- fluorine reaction. 

2. Study finite rate calculation of the same non-premixed flame. 

3. Study the interaction between mixing and chemical reaction. 

4. Study compressibility effects. 

4. Publications 

1. Hsu, A.T., “The Study of PDF Turbulence Models in Combustion,” 
9th National Aero-Space Plane Technology Symposium, November 1-2, 
1990. 

2. Hsu, A.T., “ On Recontamination and Directional Bias Problem in 
Monte Carlo Simulation of PDF Turbulence Models,” NASA CFD Con- 
ference, April 12-14, 1991, Moffett Field, California. 

3. Hsu, A.T., “Progress in the Development of PDF Turbulence Models for 
Combustion,” 10th National Aero-Space Plane Technology Symposium, 
April 23-25, 1991, Monterey, California. 

4. Hsu, A.T., “The Study of PDF Turbulence Model in Nonequilibrium 
Hydrogen Diffusion Flames” AIAA Paper 91-1780, Honolulu, Hawaii, 
June, 1991. 
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Turbulence Modeling 

Tsan-Hsing Shih 


1. Motivation and Objectives 

(1) Examine the performance of existing two-equation eddy viscosity mod- 
els and develop better models for the near-wall turbulence using direct 
numerical simulations of plane channel and boundary layer flows. 

(2) Use the asymptotic near-wall behavior of turbulence to examine the 
problems of current second-order closure models and develop new models 
with the correct near-wall behavior. 

(3) Use Rapid Distortion Theory to analytically study the effects of mean 
deformation (especially due to pure rotation) on turbulence, obtain 
analytical solutions for the spectrum tensor, Reynolds stress tensor, 
anisotropy tensor and its invariants, which can be used in the turbu- 
lence model development. 

(4) Explore the potential of the renormalization group (RNG) theory in 
turbulence modeling. 

(5) Modeling of compressible turbulent flows. 

(6) Modeling of bypass transition. 

2. Work Accomplished and Ongoing Work 
2.1 k-e model 

The k-e model is still the most widely used model for computing engi- 
neering flows. We have examined the near- wall behavior of various eddy 
viscosity models proposed by different researchers, and have studied the 
near-wall behavior of the terms in the ^-equation budget. We found that 
the modeled eddy viscosity in many existing k-e models does not possess 
correct near- wall behavior and the pressure transport term in the ^-equation 
is not modeled appropriately. Based on the near- wall asymptotic behavior of 
the eddy viscosity and the pressure transport term in the ^-equation, a new 
set of improved closure models has been obtained. In addition, a modeled 
equation for the dissipation rate is derived more rationally. This work is 
reported in NASA TM 103221 ICOMP-90-16W. 

In addition, all the existing two-equation models (except Jones & Laun- 
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der model, which unfortunately does not work well even for some simple 
flows) have an “unacceptable” wall distance parameter (y + ) in their eddy 
viscosity damping function / M (y + ). This will result in an unphysical zero 
eddy viscosity near the separation region. In addition, y + can not be well 
defined in many flows with complex geometry. To remove this deficiency, 
Dr. V. Michelassi, Dr. A. Hsu and I proposed two new eddy viscosity damp- 
ing functions, and both of them are independent of the wall distance. The 
new models have been satisfactorily tested in channel and boundary layer 
flows. This work is reported in two papers: AIAA-91-0611^ and NASA 
TM/ICOMP/CMOTTW. 

2.2 Second order modeling of near-wall turbulence 

The main emphasis is on developing a near-wall turbulence model for 
the velocity pressure gradient correlation and the dissipation tensor in the 
Reynolds-stress equation. A modeled dissipation rate equation is also derived 
more rationally. Near a wall, a reduction in velocity fluctuations normal to 
the wall becomes significant. Because of this wall effect, the viscous diffusion 
term in the Reynolds-stress equations becomes the leading term and it must 
be properly balanced by the other terms. We have used this as a model 
constraint for developing a model for the pressure and dissipation terms. To 
test the models, a fully developed channel flow and boundary layer flows 
are chosen as the test flows, for which direct numerical simulations and ex- 
perimental data axe available for comparison. The modeled Reynolds stress 
equations for the channel flow are steady one-dimensional, and for bound- 
ary layer flows are steady two-dimensional. Therefore model testing will be 
very accurate. This part of workM is reported in the paper: Proceedings 
of the International Symposium on Engineering Turbulence Modeling and 
Measurements and NASA TM 103222 ICOMP-90-0017. 

2.3 Second order modeling of a three-dimensional boundary layer 

A study of three-dimensional effects on turbulent boundary layer was 
achieved by the direct numerical simulation of a fully developed turbulent 
channel flow subjected to transverse pressure gradient (see Physics of Flu- 
ids, Vol.2 NO. 10, 1990, pp. 1846-1853). The time evolution of the flow was 
studied. The results show that, in agreement with experimental data, the 
Reynolds stresses are reduced with increasing three-dimensionality and that, 
near the wall, a lag develops between the stress and the strain rate. In ad- 
dition, we found that the turbulent kinetic energy also decreased. To model 
these three-dimensional effects on the turbulence, we have tried different 
two-equation models and second order closure models. None of the current 
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closure models can predict the reductions in the shear stress and turbulent 
kinetic energy observed in direct numerical simulations. Detailed studies 
of the Reynolds-stresses budgets were carried out. One of the preliminary 
conclusions from these budget studies is that the velocity pressure-gradient 
term in the normal stress equation ( v 2 ) plays a dominant role in the re- 
duction of shear stress and kinetic energy. These budgets have been used 
to guide the development of better models for three dimensional turbulent 
boundary layer flows. This worki 5 ! was presented in the American Physical 
Society Forty-Third Annual Meeting, November, 1990. 

2.4 The effect of rotation on turbulence 

In addition to the above studies of second order closure models, we have 
carried out some RDT analysis on simple homogeneous turbulent flows. An 
order of magnitude analysis shows that under the condition of S(q 2 )/e 
\/Rt, the equations for turbulent velocity fluctuations can be approximated 
by a linear set of equations, and if S(q 2 )/e > , then the turbulent ve- 

locity equations can be further approximated by an inviscid linear equation. 
Therefore, RDT can be used to analytically study some very basic turbulent 
flows, such as, homogeneous shear flows, irrotational strain flows and pure ro- 
tational flows. This work focuses on the effect of rapid rotation on turbulence 
using RDT. We have obtained analytical expressions for velocity, the spec- 
trum tensor, Reynolds-stress, the anisotropy tensor and its invariants. The 
solutions show that the turbulence is strongly affected by the rapid rotation. 
Using RDT, we can calculate the rapid pressure- stain term exactly and we 
can obtain very useful information for developing corresponding turbulence 
models. See the report^ for this work. 

2.5 Renormalization Group Theory (RNG) in turbulence modeling 

RN G method has been introduced to the turbulence modeling mainly in 
the Large Eddy Simulation (LES) of turbulence with a subgrid scale model. 
One also attempted to use it to develop Reynolds- averaged turbulence model 
equations, for example, k-e model equations. However, we found that there 
are a few fundamental concepts and important procedures used in the deriva- 
tion of those model equations which axe not clear and well justified. Dr. Z. 
Yang and I are working on this subject and try to explore the potential of 
RNG in the turbulence modeling. 

2.6 Modeling of compressible turbulent flows 

The turbulence models for compressible flows are of great interest in hy- 
personic flows and turbulent combustion. The modeling scheme greatly de- 
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pends on the averaging schemes (i.e., conventional average, density weighted 
average and mixed average) used in the turbulence equations. We start with 
the analysis of the turbulent equations derived from the different averaging 
schemes to see what kind of averaging scheme is most convenient for both 
turbulent modeling and applications in CFD. We concentrate on the second 
order closure model (i.e. Reynolds stress model) and two-equation model. 
Dr. W. Liou and I are working on this subject. See Reference^ for the first 
report on averaging schemes for compressible flows. 

2.7 Modeling of bypass transition 

Most common transition phenomena occurred in engineering flows axe 
bypass transition. A few papers on modeling of transition with turbulence 
models show that the bypass transition can possibly be modeled with the 
modified turbulence models developed solely for turbulent flows. However, 
most of the work in this direction was based on the parabolic two-equation 
models. We expect that the bypass transition phenomena will be more ap- 
propriately described by the elliptical equations. Then, the prediction of 
normal stresses becomes important. Because of the inability of modeling 
normal stresses with the two-equation models, we are pursuing the ellipti- 
cal Reynolds stress model equations for the bypass transition studies. Dr. 
Z. Yang and I are working on the improvement of our previous near-wall 
Reynolds stress model for the purpose of modeling bypass transition. 

2.8 Modeling of scalar turbulence: 

Modeling of scalar turbulence is of great importance in turbulent heat 
transfer. Eddy viscosity models often fail in the prediction of heat trans- 
fer in many shear flows. We have developed a set of second order closure 
models based on the joint realizability (between velocity and scalar) and the 
experiments. Dr. A. Shabbir and I are working on this subject. A paper^ 
was presented in the Lumley Symposium: Recent developments in turbulence , 
November , 1990. 


3. Publications: 

1. Shih, T.-H., 1990, “An Improved k-e Model for Near- Wall Turbulence 
and Comparison with Direct Numerical Simulation,” NASA TM 103221 
ICOMP-90-16. 

2. Shih, T.-H. and Hsu, A.T., 1991, “An Improved k-e Model for Near-Wall 
Turbulence,” AIAA-91-0611. 
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3. Michelassi, V. and Shill, T.-H., 1991, “Low Reynolds Number Two- 
Equation Modeling of Turbulent Flows,” NASA TM 104368, ICOMP- 
91-06, CMOTT-91-Ol. 

4. Shih, T.-H. and Mansour, N.N., 1990, “Modeling of Near- Wall Turbu- 
lence,” Proceeding of the International Symposium on Engineering Tur- 
bulence Modeling and Measurements , September, 1990, Dubrovnic, Yu- 
goslavia, Editors: W. Rodi, E.N. Ganic. or, NASA TM 103222 ICOMP- 
90-0017. 

5. Shih, T.-H., 1990, “Modeling of 3D Turbulent Boundary Layer Flows,” 
American Physical Society Forty-Third Annual Meeting, 1990, Cornell, 
U. Ithaca, New York. 

6. Shih, T.-H., 1991, “Rapid Distortion Theory on Homogenous Turbu- 
lence with Rapid Rotation,” CMOTT Report. 

7. Liou, W.W. and Shih, T.-H., 1991, “On the Basic Equations for the 
Second-order Modeling of Compressible Turbulence,” CMOTT-91-06. 

8. Shih, T.-H. and Shabbir, A., 1990, “Advances in Modeling the Pressure 
Correlation Terms in the Second Order Moment Equations,” the Lum- 
ley Symposium: Recent developments in turbulence, November, 1990, 
ICASE, NASA Langley Research Center, Edited by T.B. GatsJci, S.Sarkar 
and C. G. Speziale 

9. Shih, T.-H., 1990, “Advancements in Engineering Turbulence Model- 
ing,” 9th NASP Technology Symposium, Paper-105, November, 1990, 
Orlando FL. 

10. Shih, T.-H., Chen, J.-Y. and Lumley, J.L., 1991, “Second Order Mod- 
eling of Boundary Free Turbulent Shear Flows,” AIAA 91-1779. 
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Experiments and Modeling 

Aamir Shabbir 


1. Motivation and Objectives 

The usual approach in establishing the correctness and accuracy of tur- 
bulence models is to numerically solve the modeled differential equations 
and then compare the results with the experiment. However, in the case 
of a discrepancy, this procedure does not pinpoint where in the model the 
drawback lies. It is also possible that the model overcompensates one phys- 
ical phenomenon and undercompensates the other so that the net result is 
a good agreement between the two. Therefore a more desirable approach is 
to directly compare the individual terms in the equations with their mod- 
els. To achieve this objective primarily physical experiments have been used 
to carry out the second moment budgets. These can then be used to ana- 
lyze and assess various models and closure assumptions and seek improve- 
ments/modifications where models prove deficient. 


2. Work Accomplished 

2.1 Evaluation and Development of Turbulence Models for Pres- 
sure Correlations. 

A direct comparison between the pressure strain and pressure temperature- 
gradient correlations and their closure models is carried out. The flows used 
include both physical and numerical experiments on homogeneous shear flows 
and physical experiments on buoyant plumes. Models considered include 
both the linear and the more elaborate non-linear ones. It is found that 
the non-linear models provide a much better agreement with these experi- 
ments than the linear ones. A new model for the slow part of the pressure 
temperature-gradient correlation is also derived using joint realizability con- 
cept. 


2.2 On the Ratio of Mechanical to Thermal Time Scales in Tur- 
bulent Flows. 

The ratio of these time scales is very often employed in the two equation 
turbulence models. The study of Beguier et al (1978) recommended this 
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ratio to be around 2.0. The current analysis using the buoyant plume and 
homogeneous shear flow experiments shows that this value is about 3.0. It 
is shown that this departure from the commonly used value is a consequence 
of the local equilibrium assumption being not satisfied by these experiments. 

2.3 Turbulent Buoyant Transport - A Comparative Study between 
Models and Experiments. 

The more popular gradient diffusion type models for the turbulent trans- 
port (third moments) were found to underestimate the experiment by an 
order of magnitude. More complex models (Andre et al 1976, Lumley et 
al 1978), based on the simplification of exact transport equations of third 
moments, do a much better job in reproducing the third moments although 
the results are still less than satisfactory. 

2.4 Experimental Balances of Second Moment Equations for a 
Buoyant Plume. 

Despite large volume of work on second moment closure, there is very 
little experimental information available about the budgets of the second 
moments. Part of this reason stems from our inability, at present, to mea- 
sure the pressure correlations. Experimental budgets for Reynolds stresses 
and heat fluxes have been carried out for a boundary free shear flow (round 
plume) and the pressure correlations are obtained as the closing terms in 
these budgets. These budgets show how different terms in the equations are 
distributed across the flow and can be used to analyze some of the modeling 
assumptions. For example they show that the assumption of local equilib- 
rium is not justified for bulk of the flow field - an idea fundamental to the 
algebraic stress models. 

2.5 X-wire Response in Turbulent Flows of High Intensity Turbu- 
lence and Low Mean Velocities. 

This work is based on an experimental study, which was carried out at 
SUNY/Buffalo, of angular response of an x-wire, at low velocities (0.25m/s 
to lm/s). It is found that the ^-factor in the modified Cosine Law is strongly 
velocity dependent. The implications of this on multi component turbulence 
measurements are explored. Expressions are also derived for evaluating when 
the cross-flow errors begin to affect x-wire measurements. 
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2.6 Modeling of Turbomacliinery Flows using the Average Passage 
Approach. 

Turbomachineray flows are turbulent and unsteady and numerical calcu- 
lation of a flow in a multistage machine, at present, is not possible. However, 
the effects of periodic unsteadiness can be accounted through the models for 
deterministic stresses which arise in the average passage equation set ( Adam- 
czyk 1985). Exact equations governing the transport of these stresses have 
been derived and a two equation model is being developed and tested at 
present. The model uses ideas from turbulence modeling such as the gradi- 
ent diffusion type hypotheses. This work is being performed in collaboration 
with J. Adamczyk of the Lewis Research Academy, at the NASA Lewis Re- 
search center. 


3. Future Plans 

3.1 Study the effect of buoyancy on turbulence by computating flows using 
turbulence models. In addition to environmental flows, such a work also 
has industrial applications e.g. cooling of nuclear reactors and electronic 
components, and “geyser” formation in fuel tanks in microgravity. 

3.2 Seek improvements in the models for turbulent transport. In general 
the transport is not too important in most of the turbulent flows but in 
some applications, e.g. geophysical flows, the modeling of the transport 
could be critical in the success of a computation. 

3.3 Seek improvements in the existing two equation models by incorporating 
newer models for pressure correlations etc. 

3.4 Assess the models for deterministic stresses in a multistage turboma- 
chinery environment. 


4. Publications 

1. Pressure Correlations in the Reynolds Stress and Heat Flux Equations 
- A Comparison between Experiment and Models. A. Shabbir. APS 
Bulletin , Vol. 35, No. 10, Nov. 90, Abstract KC4. 

2. Evaluation of Turbulence Models for Predicting Buoyant Flows. A. 
Shabbir and D.B. Taulbee. J. Heat Transfer , 1990, Vol 112, No 4, pp 
945-951. 

3. Experiments on Round Turbulent Buoyant Plumes. A. Shabbir and 
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W.K. George. Under review for publication in J. Fluid Mechanics. 

4. Advances in Modeling Pressure Correlation Terms in the Second Mo- 
ment Equations. T.-H. Shih and A. Shabbir. Presented at the Sympo- 
sium honoring J. Lumley’s 60th birthday , November 90, NASA Langley 
Research Center. 

5. X-wire Response in Turbulent Flows of High Intensity Turbulence and 
Low Mean Velocities. A. Shabbir, P.D. Beuther, and W.K. George. 
Submitted to Experimental Thermal and Fluid Science. 


5. References 

Adamczyk, J. J. “Model Equation for Simulating Flows in Multistage 
Turbomachinery”, ASME Paper No. 85-GT-226. (1985) 

Andre J.C., G. De Moor, P. Lacarrere and R. du Vachat “Turbulence 
Approximation for Inhomogeneous Flows: Part I. The Clipping Approx- 
imation”, J. Atmos. Sci., Vol. 33, pp. 476-481, (1976) 

Beguier, C., I. Dekeyser and B. E. Launder “Ratio of Scalar and Velocity 
Dissipation Time Scales in Shear Flow Turbulence” ,” Phys. Fluids , Vol. 
21, pp. 307-310 (1978). 

Lumley, J. L., 0. Zeman and J. Siess “The influence of Buoyancy on 
Turbulent Transport”, J. Fluid Mech ., Vol. 84, pp. 581-597 (1978). 
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Modeling of Compressible Turbulent Shear Flows 

William W. Liou 


1. Motivation and Objectives 

Despite all the recent development in computer technologies and numer- 
ical algorithms, full numerical simulations of turbulent flows axe feasible only 
at moderate Reynolds numbers and for flows with relative simple geometries. 
Turbulence models provide alternatives in the pressing need for the prediction 
of turbulent flows and, in fact, have become an important pacing factor for 
the successful development of computational fluid dynamics (CFD). With 
the advent of supercomputers, however, it has become more affordable to 
apply second order closure models in the prediction of flows with complex 
effects; such as strong curvature, three-dimensionality and compressibility. 
The main goal of this research is to develop new second order moment clo- 
sures for compressible turbulence. It has been shown that the models based 
on the extension of those developed originally for incompressible flows fail 
to predict adequately turbulent flows at high Mach numbers. In this at- 
tempt, the compressibility effects will be explicitly considered. A successful 
development of these models that take into account directly the compress- 
ibility effects may have a range of technological implications in the design of 
supersonic and hypersonic vehicles. 


2. Work Accomplished 

During this early stage of the task, the goal is to obtain an objective yet 
comprehensive understanding of the development and the current status of 
compressible turbulence modeling. Due to the variable density effects in com- 
pressible flows, density correlation terms appear in the governing equations 
for the mean flow, if the conventional ensemble averaging technique is ap- 
plied. These terms do not exist in incompressible flows and need to be mod- 
eled. On the other hand, the mass- weighted- averaging or Favre procedure 
generates a set of mean equations that have the similar forms as they are for 
incompressible flows. One may then incline to use the incompressible analog 
in compressible flow calculations. A simple test, however, should show that a 
direct application of the exact incompressible models fails. One of the main 
effects that is excluded in incompressible models is the finite propagation 
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speed of disturbances. In compressible flows, modulation of flow properties 
occurs only within Mach cones of influence, with acoustic time delay. This 
introduces additional scales for the transport properties. Caution also needs 
to be exercised in comparing Favre-averaged calculations with experiments, 
since the differences between Favre-averaged and measured quantities may 
not be negligible at high Mach numbers. 

Compressible turbulence modeling is still in its infancy. This appears 
to be true both theoretically and computationally. Recently, a new concept 
called dilatation dissipation was proposed. Dilatation dissipation, as op- 
posed to the solenoidal dissipation in incompressible flows, accounts for the 
viscous dissipation of turbulent kinetic energy due to volume fluctuations. 
Models accommodating this effect show the importance of this additional 
drain of turbulent kinetic energy in order to obtain adequate model predic- 
tions. Dilatation dissipation appears to be among the direct consequencies 
of compressibility effects. 

Another school of thought on compressibility effects focuses on the 
changes of turbulence structures at high Mach numbers. Note that these 
structures are identified by using conditional sampling techniques in exper- 
iments. Due to the communicability problem between interacting elements, 
structures that are highly efficient in extracting energy from the mean flow at 
low Mach numbers no longer prevail as the Mach number increases. They are 
replaced by structures that are less sensitive to thq Mach number. This se- 
lective amplifying behavior is describable by quasi-linear theory, which view 
the turbulence energetics as physical manifestations of ongoing nonlinear 
instability in turbulent shear flows. 

The above mentioned matters are described in detail in an ICOMP / 
CMOTT report [1] that is in preparation. Equations for the second order 
moments and the mean flow as a result of the application of different averages 
will also be given. Modeling methodologies used in compressible flow calcula- 
tions will be reviewed. The evaluation process performed during the present 
stage of the research has identified avenues that will be pursued during the 
next period of this task. 


3. Future Plans 

(1) Develop second order models for compressible turbulence based on en- 
semble averages. This may be assisted by first developing k — e types of 
models to identify important mechanisms. 

(2) Develop unconventional models that incorporate explicitly the charac- 
teristics of the structures of compressible turbulence. 


i 
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The developed models will be applied to certain benchmark flows with and/or 
without chemical reactions. 


4. Publications 

1. Liou, W. W. and Shih, T.-H., “On the Basic Equations for the Second- 
order Modeling of Compressible Turbulence,” CMOTT-91-06, 1991. 

2. Liou, W. W. and Morris, P. J., “An Comparison of Numerical Methods 
for the Rayleigh Equation in Unbounded Domains,” CMOTT-91-05, 
1991. 
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RNG in Turbulence and Modeling of Bypass Transition 

Zhigang Yang 

1. Motivation and Objectives 

Since I joined CMOTT on July 1990, I have been working on two 
research projects. The first project concerns the Renormalization Group 
(RNG) analysis of turbulence and the second project is on the calculation of 
bypass transition through turbulence modeling. In addition, the preparation 
of two papers on work performed at was completed. 

Application of RNG in turbulence was proposed by Yakhot and Orszag 
in 1986. RNG is a process which eliminates systematically the small scales, 
and represents the effect of those eliminated small scales on the uneliminated 
large scales as the changes in the transport properties. It is because of this 
property of RNG that Yakhot and Orszag suggested that RNG could be 
used as a model builder in turbulence modeling. They also presented a 
k — e model in their 1986 paper. However, this paper is lengthy, with many 
unstated assumptions. Our aim is first to understand, and to validate the 
RNG approach in turbulence through an independent study. We will then 
study the possibility of constructing RNG based turbulence models, and try 
to proceed to do the turbulence modeling through RNG in parallel with the 
classical approach. We will also compare the numerical predictions made by 
RNG models and by classical models against data from Direct Numerical 
Simulation and against experimental data from different benchmark cases. 

In a quiescent environment, the transition is initiated by the instability 
of the laminar boundary layer to Tollmien-Schlichting waves. These waves 
are amplified with streamwise distance and eventually breakdown into tur- 
bulent spots, which are precursors of turbulent boundary layer. While in an 
environment with high freestream turbulence, the transition is found to be a 
bypass one in which turbulent spots are formed without Tollmien-Schlichting 
wave amplification. The formation of turbulent spot is a random process, 
and flow within a turbulent spot is almost fully turbulent. This suggests the 
possibility of using turbulence modeling to describe and predict the bypass 
transition. There have been some works in this direction, primarily using 
different versions of two equation models. Bypass transition is predicted, as 
the level of the freestream turbulence is increased. However, it is found that 
the predicted transition is much sharper than that observed in the experi- 
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ment. In addition, the predicted transition depends on the description of the 
initial profiles to an certain extent. The works we propose to do are twofold. 

1) We will be using a low Reynolds number version of Reynolds stress model 
rather than the two equation model. This would bring in more physics, and 
hopefully would fit better the complicated flows such as bypass transition. 

2) We will be using an elliptic solver rather than a parabolic solver for this 
boundary layer transition. This way, we will be able to include the effect of 
the leading edge. The testing case will be flow passing a flat plate. Both the 
zero pressure gradient case and the non- zero pressure gradient case will be 
tested. 

2. Work Accomplished 

1. Nonlinear dynamics near the stability margin in rotating pipe flow. 

The nonlinear evolution of marginally unstable wave packets are studied 
in rotating pipe flow. These flows depend on two control parameters, which 
may be taken to be the axial Reynolds number Re and the rotation rate q. 
Marginal stability is realized on a curve in the (Re, q ) plane, and we explore 
the entire marginal stability boundary. As the flow passes through any point 
on the marginal stability curve, it undergoes a supercritical Hopf bifurcation 
and the steady base flow is replaced by a traveling wave. The envelope of 
the wave system is governed by a complex Ginzburg- Landau equation. The 
Ginzburg-Landau equation admits Stokes waves, which correspond to stand- 
ing modulations of the linear traveling wavetrain, as well as traveling wave 
modulations of the linear wavetrain. Bands of wavenumbers are identified in 
which the nonlinear modulated waves are subject to a side-band instability. 

This work[l] was reported in APS/DFD meeting in November 1990. A 
paper for this work has been submitted to JFM for consideration of publica- 
tion. The paper is co-authored with S. Leibovich of Cornell University. This 
work was supported by AFOSR-89-0346. 

2. Unstable viscous wall modes in rotating pipe flow. (AIAA Paper No. 
91-1801) 

Linear stability of flow in rotating pipe is studied. These flows depends 
on two parameters, which can be taken as the axial Reynolds number Re and 
the rotating rate q. In the region of Re > 1 and q = 0(1), the most unstable 
modes are wall modes. The wall modes are found to satisfy a simpler set of 
equations containing two parameters rather than four parameters as in the 
full linear stability problem. The set of equations is solved numerically and 
asymptotically over a wide range of the parameters. In the limit of Re — > oo, 
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the eigenvalue reaches the inviscid limit and the eigenfunction shows a two 
layer structure. The eigenfunction reaches the inviscid limit over the main 
part of the domain, while near the wall of the pipe, the eigenfunction is 
represented by a viscous solution of the boundary layer type. 

This work[2] is to be presented at the AiAA 22nd Fluid Dynamics, 
Plasma Physics and Lasers Conferences, June 24-26, 1991. The paper is co- 
authored with S. Leibovich of Cornell University. This work was supported 
by AFOSR-89-0346. 

3. RNG in turbulence modeling. 

This work is done in collaboration with Dr. T.H. Shih of CMOTT. In 
this work, we carry out an independent study of the work done in the paper 
by Yakhot and Orszag up to their derivation of k — e equation. Many of their 
results are repeated. However, we also found some discrepances, and some 
unclaimed assumptions in their derivations. 


3. Future Plans 

1. Currently, we are testing and improving the low Reynolds number 
version of Reynolds stress model proposed by Shih and Mansour (1990) in 
the simple shear flows, such as channel flow and boundary layer flow. This 
model is going to be used in the calculation of bypass transition. 

2. We will carry an independent derivation of k — e equation using RNG, 
and compare it with the one presented by Yakhot and Orszag. We will also 
compare the prediction of RNG k — e model with the other k — e models, in 
both the high Reynolds number case and the low Reynolds number case. 


4. Publications 

[1] Yang, Z. and Leibovich, S. 1990 “Nonlinear dynamics near the stability 
margin in rotating pipe flow”, Submitted to J. Fluid Mech. 

[2] Yang, Z. and Leibovich, S. 1990 “Unstable viscous wall modes in rotating 
pipe flow”, AIAA Paper 91-1801. 

5. References 

Yakhot, V. and Orszag, S.A. 1986 Renormalization group analysis of 
turbulence. J. Sci. Comput. Vol. 1, No. 1, 3-51. 

Shih, T.H. and Mansour, N.N. 1990 Modeling of near wall turbulence. 
NASA TM-103222, ICOMP-90-0017. 
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Turbulence Modeling in Supersonic Combusting Flows 


Tawit Chitsomboon 


1. Motivation and Objectives 

To support the National Aerospace Plane project, the RPLUS3D CFD 
code has been developed at the NASA Lewis Research Center. The code 
has the capability to solve three dimensional flowfields with finite rate com- 
bustion of hydrogen and air. The combustion processes of the hydrogen- 
air system are simulated by an 18-reaction path, 8-species chemical kinetic 
mechanism. The code uses a Lower-Upper (LU) decomposition numerical 
algorithm as its basis, making it a very efficient and robust code. Except 
for the Jacobian matrix for the implicit chemistry source terms, there is no 
inversion of a matrix even though it uses a fully implicit numerical algorithm. 

The main purpose of this work is to incorporate a k-e (two equation) 
turbulence model into the RPLUS3D code. 


2. Work Accomplished 

Since February 1990, when this work was started, some of the more 
important accomplishments axe categorized as follows: 

1) Add a k-e turbulence model: The model selected is in a high Reynolds 
number form. The low Reynolds number form could not be used economi- 
cally in the case of a three dimensional flow with chemical reaction since it 
demands too much in computer resources. The addition was designed to be 
as modular as possible but some interactions with the main code are needed 
in order to be more efficient. The first test case tried was a Mach 0.5 flow 
over a flat plate. The velocity profile compare very well with the log-law 
profile. The friction coefficient also compares well with the Van Driest corre- 
lation. More validations will be performed for other flows such as free shear 
layer and jet flows. 

2) Improved accuracy and convergence rate: According to a stability 
analysis of a model equation, it is shown that the RPLUS3D code exces- 
sively added artificial damping to the right hand side of the algorithm while 
at the same time it overestimated the spectral radii on the left hand side. 
These excessive additions would not give an optimum convergent rate. Mod- 
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ifications were made such that true directional spectral radii were added to 
the left hand side and the artificial damping terms were reduced to optimum 
values in accord with the stability analysis. 

A test run was made of a Mach 4 flow of air over a 10 degree compression 
ramp. It was found that the modified code converged to machine zero about 
five times faster than the original code while at the same time it somewhat 
improved the shock resolution. 

3) Added consistent damping terms at block interfaces: It had been ob- 
served that at the block interface of a multi-block grid, wiggles developed. 
This happened because the damping terms at the interface were not consis- 
tent with those at the interior points. 

4) Validated RPLUS3D with laminar flows over flat plates: The test 
cases run were for Mach numbers of 0.1, 0.3 and 0.5. Results of all cases 
seemed to be good except for the region of high curvature of the velocity 
profile near the edge of the boundary layer. 

5) Add two dimensional capability to the code: It turned out that this 
is not a trivial task especially for a finite volume code like RPLUS3D. It is 
nice that now the code can solve a 2D flow without having to carry the 3rd 
direction along as a redundancy. There is, of course, no need to maintain a 
separate 2D code. 

6) Implemented a local time stepping capability: The original code al- 
ways ran at a fixed time step of 1 second. For most flows this corresponds to 
using a very large CFL number which may not be conducive to a fast con- 
vergence rate. With the local time stepping, it was found that an optimum 
CFL number was, in agreement with other investigators, around 5 to 7. 

7) Implemented implicit boundary conditions: This addition enhanced 
the convergence rate by about 30 percent at the expense of a more complex 
code and an increase of about 20 percent in CPU time per iteration step. 
There seems to be, then , no net advantage of the implicit boundary condition 
except maybe in the area of robustness. 

8) Changed input file: Instead of having to scan a whole subroutine to 
set up a problem, a user now can change numbers in a small file of length 
about one page. To start up a run from a previous run one now needs only to 
change a parameter in this input file without having to recompile the input 
subroutine as before. 

9) Solved two 3D hypermixing flow fields of W.Hingst and D. Davis: This 
work was performed in collaboration with Dr. A.C. Taylor of Old Dominion 
University. My task was to set up the program for these particular problems 
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and to generate the grid for one of the problems. 


3. Future Plans 

The work plan for the year 1991 consists of both basic model implemen- 
tations and practical applications of the code : 

• Continue to validate the baseline k-e model. 

• Add compressibility effects to the base line model. 

• Apply the base line model to re-solve the 3D flowfields mentioned in the 
previous section. 
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Development of New Flux Splitting Schemes 

Meng-S. Liou and Christopher J. Steffen, Jr. 


1. Motivation and Objectives 

Maximizing both accuracy and efficiency has been the primary objective 
in designing a numerical algorithm for computational fluid dynamics (CFD). 
This is especially important for solution of complex 3D systems of Navier- 
Stokes equations which often include turbulence modeling and chemistry 
effects. Recently, upwind schemes have been well received for both their 
capability of resolving discontinuities and their sound theoretical basis in 
characteristic theory for hyperbolic systems. With this in mind, we present 
two new flux splitting techniques for upwind differencing. 


2. Work Accomplished 

The first method is based upon High-Order Polynomial Expansions 
(HOPE) of the mass flux vector [1]. The present splitting results in positive 
and negative mass flux components that vanish at M—0. Thus the error in 
the Van Leer scheme which results in the diffusion of the boundary layer is 
eliminated. We also introduce several choices for splitting the pressure and 
examine their effects on the solution. 

The second new flux splitting is based on the Advection Upwind Split- 
ting Method (AUSM for short) [2]. In Navier-Stokes calculations, the diffu- 
sion error present in Van Leer’s flux splitting scheme corrupts the velocity 
vector near the wall. In the AUSM, a proper splitting of the advective ve- 
locity component leads to an accurate resolution of the interface fluxes. The 
interface velocity is defined using the Mach number polynomial expansion 
in the mass flux, then the convective fluxes follow directly. Again, several 
choices of pressure splitting are possible among which a simple Mach num- 
ber splitting according to characteristics appears to be the best in terms of 
accuracy. The scheme has yielded results whose accuracy rivals, and in some 
cases surpasses that of Roe’s method, at reduced complexity and computa- 
tional effort. The calculation of the hypersonic conical flow demonstrates 
the accuracy of the splittings in resolving the flow in the presence of strong 
gradients. The second series of tests involving the 2D inviscid flow over a 
NACA 0012 airfoil demonstrate the ability of the AUSM to resolve the shock 
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discontinuity at transonic speed and the level of entropy generation at the 
stagnation point. 

In the third case we calculate a series of supersonic flows over a cir- 
cular cylinder. The Roe splitting in all conditions and grids tested yielded 
anomalous solutions (sometimes referred to as the carbuncle phenomenon), 
which could appear as non-symmetric, protuberant, or indented contours. 
The AUSM gave expected solutions in all calculations. 

The fourth test deals with a 2D shock wave/boundary layer interaction. 
This provides an opportunity to accurately resolve a laminar separation re- 
gion and to compare the ability to resolve a non grid-aligned shock with 
other methods. 

3. Future Plans 

Future plans axe primarily concerned with the AUSM. A detailed sta- 
bility analysis for this new technique will be useful. The idea of splitting the 
advective velocity opens up a whole family of potential schemes. Therefore, a 
comprehensive study of the interaction between various pressure and advec- 
tion velocity splitting methods is necessary to optimize both accuracy and 
efficiency. Additionally, a 2D turbulent calculation would be a good test of 
the scheme’s ability to solve a coupled system of k — e equations. 


4. Publications 

1. Liou, M.-S. and Steffen, C.J.Jr., “High-Order Polynomial Expansions 
(HOPE) for Flux- Vector Splitting,” (to be presented at the ICES’91 
Conference, August 11-16, 1991) 

2. Liou, M.-S. and Steffen, C.J.Jr., “Development of a New Flux Splitting 
Scheme,” (to be presented at the AIAA Tenth CFD Conference, June 
24-26, 1991) 
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Progress of Simulations for Reacting Shear Layers 

Sheng-Tao Yu 


1. Motivation and Objectives 

In the past six months, the effort was devoted to the development of 
a high speed, chemically reactive shear layer test rig. The purpose of the 
experiment is to study the mixing of oxidizer and fuel streams in reacting 
shear layers for various density, velocity, and Mach number. The primary 
goal is to understand the effects of the compressibility upon mixing and 
combustion in a fundamental way. Therefore, a two-dimensional shear layer 
is highly desirable for its simplicity to quantify the compressibility effects. 

The facility consists of a two-stream wind tunnel with two indepen- 
dent gas supplies. After passing through flow-management devices located 
upstream, each gas stream expands to its predetermined Mach number by 
means of a contoured center body and tunnel walls. Various combinations 
of flow conditions of high-speed stream and low-speed stream allows for the 
systematic study of mixing and reactions of compressible shear layers. 


2. Work Accomplished 

The RPLUS 2D code is used to calculate the flow fields of different 
sections of the test rig. The emphasis was on the supersonic nozzle design, 
the vitiation process for the hot air stream and the overall thermodynamic 
conditions of the test matrix. 

The k — e turbulence model with wall function has been successfully 
implemented in the RP LU S code. The k and e equations are solved simulta- 
neously and the the LU scheme is used to make it compatible with the flow 
solver. The coupling between the flow solver and the k — e solver depends 
on the turbulence viscosity only, and the k — e solver is separated from the 
flow solver to reduce the complexity. The newly developed code has been 
used for the compressible shear layer calculations. Many cases of the com- 
pressible free shear layer with various convective Mach numbers and density 
ratios have been simulated using the compressible k — e solver. The results 
are summarized in two technical papers. 5,6 Currently, the k — e solver is a 
standard feature in the RPLUS 2D code and the code has been distributed 
to the industry and universities through NASP group. Locally, Duncan and 
Tsai are using the k — e solver for their research work. 
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3. Future Plans 

Physical phenomena of the reactive free shear layer can not be ade- 
quately described by the k — e model coupled with the Reynolds averaged 
flow equations. The properties of the vortical flows which dominate the 
whole flow field of free shear layers can only be illustrated by time marching 
numerical method with accurate spatial resolution. Traditionally, the spec- 
trum methods were used for this kind of applications. However, very limited 
success has been reported for compressible flows using spectrum methods. 
On the other hand, recent development show promising results using high 
order central differencing and Essentially Non-Oscillatory (ENO) schemes. 
One developed by Lele of Stanford university using high order compact dif- 
ferencing is especially interesting. Future work includes Direct Numerical 
Simulation (DNS) of Navier Stokes equations for the chemically reactive 
flow and application to free shear layers. 

In additional to the above mentioned work, I will serve as a consultant 
for the k — e solver in the RPLUS code. 

4. Publications 

1. S. T. Yu, J. S. Shuen, and Y-L P. Tsai, “Three Dimensional Calculations 
of Supersonic Combustion Using a LU Scheme,” to appear in the J. of 
Comput. Phys. 

2. S.T. Yu, C.L. Chang, and C.L. Merkle, “Solar Rocket Plume/Mirror 
Interactions, ” submitted to the J. of Spacecraft and Rockets. 

3. S.T. Yu, B.J. McBride, K.C. Hsieh, and J.S. Shuen, “Hypersonic Flows 
Simulations with Equilibrium or Finite Rate Chemistry,” submitted to 
Computers & Fluids for publication 

4. S.T. Yu, “A Convenient Way to Convert 2D CFD Codes to Axisym- 
metric Ones,” submitted to J. of Propulsion and Power as a technical 
note. 

5. S.T. Yu, C.D. Chang, and C.J. Marek, “Modem CFD applications for 
the Design of a Reacting Shear Layer Facility,” presented at the AIAA 
Science Meeting, 1991. 

6. S.T. Yu, C.D. Chang, and C.J. Marek, “Simulation of Free Shear Layers 
Using a Compressible k-e Model,” accepted for presentation at the AIAA 
Propulsion Conference, 1991. 
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Dr. Meng-Sing Liou 
Senior Scientist 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 

Technical Leader 

Dr. Tsan-Hsing Shih 

Senior Research Associate 

ICOMP, NASA Lewis Research Center 
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Dr. Marvin E. Goldstein 
Chief Scientist 

NASA Lewis Research Center 

Professor John L. Lumley 
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Aerospace Engineering 
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Professor Eli Reshotko 
Department of Mechanical and 
Aerospace Engineering 
Case Western Reserve University 
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Member Listing 


Names/Term 

Affiliation 

Research Areas 

Chitsomboon, Tawit 
5/90 - present 

ICOMP 

Algorithms, 

Code Development 

Hsu, Andrew T. 
5/90 - present 

Sverdrup Tech., Inc. 

PDF Turbulence Modeling 

Lang, Nancy J. 
5/90 - 8/90 

NASA Intern 
Univ. of Michigan 

Turbulence Modeling, 
CFD 

Liou, Meng-Sing 
5/90 - present 

NASA LeRc 

CFD Algorithms, 
High-speed Flow 

Liou, William W. 
11/90 - present 

ICOMP 

Compressible Flow 
Modeling, Weakly 
Nonlinear Wave Models 

Shabbir, Aamir 
5/90 - present 

NRC 

Buoyancy Effects on 
Turbulence 

Shili, Tsan-Hsing 
5/90 - present 

ICOMP 

Turbulence Modeling 

Steffen, Christopher J. Jr. 
10/90 - present 

NASA LeRc 

Upwind Algorithms 
Incompressible Flow 

Yang, Zhigang 
7/90 - present 

ICOMP 

RNG in Turbulence 
Modeling, Modeling of 
Bypass Transition 

Yu, Sheng-Tao 
3/90 - present 

Sverdrup Tech., Inc. 

Modeling of Chemical 
Reacting Flows, 

Direct Numerical Simulations 
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Appendix B 

Seminars and Technical Meetings 


B.l CMOTT Seminar Seminars 

The purpose of these seminars is to exchange ideas and opinions about 
the latest developments and current state of turbulence and transition re- 
search. The speakers are invited from within and outside of the NASA 
LeRC, including foreign speakers. This seminar series complements the in- 
formal CMOTT technical group meetings. 

The abstracts of the seminars are given below. 
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INSTITUTE FOR COMPUTATIONAL 
MECHANICS IN PROPULSION 


ICOMP SEMINAR SERIES 


Measurements in a Circular Jet of Helium into 
Quiescent Air* 

by 

John L. Lumley 
Cornell University 

Monday, June 4, 1990 
9:00-10:30 a.m. 

ERB Building 5, Room 119 


The ability to model an inert flow with density fluctuations allows one to model a 
diffusion flame with fast chemistry, and gives insight into one aspect of the 
dynamics of compressible flows. With this in view, these measurements were 
made to calibrate a second order model for such flows. Measurements using 
shuttle-mounted hotwire and Way-Libby probes are described. For comparison, 
a jet of air into air was also measured. Moments up to fourth order were 
computed. The data are used to evaluate various modeling assumptions. 
Attempts are made to explain the greater spreading rate of the helium jet. 

‘From the Ph.D. thesis of N. R. Panchapakesan 


Contact: Charles Feiler, PABX 3-6681 


30 



CENTER FOR MODELING OF 
TURBULENCE AND TRANSITION 

CMOTT SEMINAR SERIES 


ENGINEERING TURBULENCE MODELING 
Present and Future? 

Tsan-Hsing Shih 

Center for Modeling of Turbulence and Transition 


Wednesday, July 11, 1990 
3.00-4:00 PM 

ERB Building 5, Room 119 


A summary of the present position of eddy-viscosity models (e.g. 
k — e) and second-order closure models (Reynolds stress models) 
is presented. Typical examples (comparisons between model predic- 
tions and experiments) show their abilities as well as their limitations. 

Development of more advanced and complex schemes is discussed. 
The inclusion of such models into a CFD commercial code is feasible, 
but need intensive work - a cooperative effort! 


Contact: T.-H. Shih, PABX 3-6680 
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TURBULENCE AND TRANSITION 
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FUNDAMENTAL IN-HOUSE EXPERIMENTS 
TO SUPPORT THE LEWIS TURBULENCE 
MODELING PROGRAM 

Edward J. Rice 

Inlet, Duct and Nozzle Flow Physics Branch 

Wednesday, July 25, 1990 
3.00 PM 

Building 5, Room 119 

The fundamental experiments conducted by members of the Inlet, Duct and Nozzle 
Flow Physics Branch of IFMD have been primarily for support of the shear flow control 
effort within the Division. Although this support effort is expected to continue, some work 
could be diverted (and in time expanded) to support the Turbulence Modeling Program 
currently underway within ICOMP. This new emphasis can be accomplished through an 
interactive effort between the Numerical Analysts and the Experimentalists. The current 
and planned experiments will be summarized in this talk and include: unsteady flow around 
airfoils (stationary and oscillating), 2D rapid diffusers (backward facing ramp), aspirated 
backward facing step, circular and rectangular jets (subsonic and supersonic), duad stream 
supersonic shear layer (annular geometry), and boundary layer transition. A swirl generator 
within the plenum of the CW-17, ERB rig allows the addition of swirling flow to any of 
the jet or shear layer experiments. The available and planned experimental instrumentation 
include: single, X-wire and three-wire hot wire anemometry, single and two element corona 
probe, multiple microphone and pressure transducer channels, Schlieren and laser-sheet flow 
visualization, and conventional and fiber-optic two-component LDA systems. The new 16 
channel anemometry system allows simultaneous measurement with 8 X-wires. Some sample 
data will be presented to illustrate the current capability. 

This is intended to be a group discussion type of workshop. The presentation will take 
30 minutes and will be followed by 30 minutes of discussion between the numerical and 
experimental participants. The presentation will be general in nature intended mainly to 
acquaint the numerical analyst with the experimental capability which can support the 
numerical program through a cooperative effort. 

For additional information contact: T.-H. Shih, PABX 3-6680 
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CENTER FOR MODELING OF 
TURBULENCE AND TRANSITION 


CMOTT SEMINAR SERIES 


CFD-RELATED RESEARCH AT THE 
UNIVERSITY OF BRUSSEL 
AND TURBULENCE RESEARCH ACTIVITIES IN 

EUROPE 


Prof. Charles Hirsch 
Vrije Universiteit Brussel 


Wednesday, July 25, 1990 
1.00 - 2.00 PM 
Building 5, Room 119 


For additional information contact: L. Povinelli, PABX 3-5818 
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TURBULENCE AND TRANSITION 

CMOTT SEMINAR SERIES 


CHAOTIC WANDERINGS THROUGH A LAND 
OF TURBULENCE MODELS 


Russ Claus 

Aerothermocliemistry Branch 


Wednesday, August 8, 1990 
3.00 PM 

Building 5, Room 119 

A turbulence model user details his pursuit of the mythological “correct 
answer”. The talk opens with some failed attempts to calculate seperated 
flows using a two equation turbulence model. Following this, additional 
attempts to achieve the “correct answer” through a series of Direct Nu- 
merical Simulations and Large Eddy Simulations will be described and the 
limitations of these approaches axe highlighted. Finally, a re-examination 
of two-equation and a second order closure calculations of a jet in crossflow 
will be discussed. 

A brief discussion of some turbulence modeling efforts being supported 
under the NASA SBIR program will also be described. This includes RNG 
modeling with Orszag and Yakhot, and PDF modeling for compressible 
flows with Kollman and Farshchi. 

For additional information contact: T.-H. Shih, PABX 3-6680 
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TURBULENCE AND TRANSITION 

CMOTT SEMINAR SERIES 


A COMPARATIVE ANALYSIS OF TWO 
EQUATION TURBULENCE MODELS 

Nancy Lang 
CFD Branch 


Wednesday, August 22, 1D90 
3.00 PM 

Building 5, Room 119 

Several two equation models have been proposed and tested against 
benchmark flows by various researchers. For each study, different numerical 
methods or codes were used to obtain the results which were an improve- 
ment or success over some other model. However, these comparisons may 
be overshadowed by the different numerical schemes used to obtain the re- 
sults. With this in mind, several existing two equation turbulence models, 
including k-e and k-r models, are implemented into a common flow solver 
code for near wail turbulent flows. Calculations are carried out for low 
Reynolds number, two dimensional, fully developed channel and boundary 
layer flows. The accuracy of the different models is established by com- 
paring the turbulent kinetic energy, mean velocity, and shear stress profiles 
with the direct numerical simulations and experimental data. 

For additional information contact: T.-H. Shih, PABX 3-66S0 
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ON LENGTH SCALE EQUATIONS IN 
TWO-EQUATION TURBULENCE MODELS 

Micha Wolfshtein 
Faculty of Aerospace Engineering 
Technion, Israel Institute of Technology 
Haifa, Israel. 


Wednesday, August 29, 1990 
3.00 PM 

Building 5, Room 119 

Two-equation turbulence models show a great survivability between the very successful 
mixing length models and the Reynolds stress models. Yet the models are often critisized, 
mainly on the wide scatter of the computed results and on the difficulties encountered in the 
derivation of a reliable scale equation. The problem has been difficult to resolve due to: 

(i) The large resources required for developing solvers and to run test cases on computers; 

(ii) The theoretical difficulties to derive turbulence models from the Navier Stokes equations. 
The demands from a ”good n turbulence model are difficult to satisfy, and are often conflicting 
with one another. This point will be illustrated in a discussion of some possible approaches 
to this problem. In particular we shall refer to well established models like the dissipation 
or length scale models, as well as newer models like the volume of turbulence or time scale 
models. A generalised two-equation turbulence model will be used to demonstrate a possible 
approach for the improvement of two-equation models. Finally, a fourth order boundary layer 
solver for the generalized two equations model will be presented. The solver can handle both 
compressible and incompressible flows, with any two-equation model, with or without wall 
functions. Some results will be presented, for a flat plate boundary layer and for unseparated 
diffuser flows. 


For additional information contact: T.-H. Shih, PABX 3-66S0 
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EFFECT OF ACOUSTIC EXCITATION ON 
STALLED FLOWS OVER AN AIRFOIL 

Khairul Zaman 

Inlet, Duct and Nozzle Flow Physics Branch 


Wednesday, September 5, 1990 
3.00 PM 

Building 5, Room 119 

Experimental results on the subject are to be s umm arized focussing atten- 
tion on post-stalled flows, i.e. flows that are fully separated from near the 
leading edge of the airfoil. The excitation results in a tendency towards 
reattachement, which is accompanied by an improved airfoil performance, 
although the flow may still remain fully seperated. It is observed that with 
increasing excitation amplitude, the effect becomes more pronounced but 
shifts to a Strouhal number which is much lower than that expected from 
linear, inviscid instability of the seperated shear layer. In addition, some 
results from a recent experiment on supersonic jets will be briefly reviewed 
emphasizing need for collaboration between experiment and computations. 


For additional information contact: T.-H. Shih, PABX 3-6680 
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COMPUTATION OF ELLIPTIC FLOWS USING 
LOW-REYNOLDS-NUMBER TWO-EQUATION 

MODELS 

Vittorio Michelassi 
Energy Engineering Department 
University of Florence, Italy 

Wednesday, September 19, 1990 
3.00 PM 

Building 5, Room 119 

Two new low-Reynolds number forms of the k — e model will be presented. 
These exhibit better stability and stiffness characteristics as compared to 
the previous formulations. Model generality is improved by formulating 
the damping functions so that they do not depend on the wall distance. 
The proposed formulations are compared with eight other low Reynolds 
number two-equation turbulence models by computing the fully developed 
channel flow and the incompressible flow past a hill. Results are compared 
with available direct numerical simulation and experimental data. The flow 
solver is based on the approximate factorization technique and the artificial 
compressibility method requiring no (or very little) numerical damping. A 
simple linearization technique for the turbulence model source terms based 
on Taylor series expansion ensures implicit algorithm stability for all the 
models tested. Both numerical accuracy and computational efficiency are 
discussed. 



For additional information contact: T.-H. Shih, PABX 3-6680 
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Evaluation of Models for the Pressure Terms 
in the Second Moment Equations 

Aamir Shabbir 

Wednesday, October 10, 1990 
3.00 PM 

Building 5, Room 119 

The usual approach in establishing the correctness and accuracy of turbur 
lence models is to numerically solve the modeled differential equations and 
then compare the results with the experiment. However, in the case of a 
discrepancy, this procedure does not pinpoint where the drawback of the 
model lies. It is also possible that the model overcompensates one physical 
phenomenon and undercompensates another resulting in good agreement 
between the prediction and experiment. A more desirable approach is to 
individually compare each term in the equations with its model. This talk 
will focus on such a comparison for the pressure correlations appearing in 
the second moment equations. At present these correlations can not be mea- 
sured but can be obtained by balancing the turbulent flux and Reynolds 
stress equations. Using this approach, pressure correlations will be directly 
compared with their models, which range from simple linear to more elab- 
orate nonlinear ones. Also results will be shown for a newly developed 
model for the return to isotropy part of the pressure temperature-gradient 
correlation. The data used is from the homogeneous shear flow and the 
buoyant plume experiments as well as from the direct numerical simulation 
of homogeneous shear flow. 

For additional information contact: T.-H. Shih, PABX 3-6680 
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HEAT TRANSFER WITH INHOMOGENEOUS 
FREE STREAM TURBULENCE 

S.N.B. Murthy 

School of Mechanical Engineering 
Purdue University, West Lafayette, IN 47907 


Thursday, October 25, 1990 
2:00 PM 

Building 5, Room 119 

The presence of free stream turbulence (FST) in a wall-bounded flow with 
heat transfer presents several interaction complexities. These depend upon 

(a) the state of the boundary layer, laminar, fully turbulent or transitional, 

(b) the nature of FST, especially its “peakiness” or inhomogeneity, and (c) 
other complications such as pressure gradient, geometry and cooling. An 
investigation is in progress on the possible application of (a) large eddy in- 
teraction hypothesis (based on Lumley’s rational description of turbulence) 
and (b) spectral analogy between heat and turbulence kinetic energy to the 
simplest case of a flat wall, fully developed turbulent boundary layer with 
zero pressure gradient when there is heat transfer in the presence of homo- 
geneous FST. An extension of the approach to the case of other types of 
boundary layers with inhomogeneous FST will be discussed. 


For additional information contact: T.-H. Shih, PABX 3-66S0 
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PRESSURE-STRAIN MODELLING AND 
THE RETURN TO ISOTROPY 


M.M. Gibson 

Imperial College of Science Technology, 
London, England. 


Wednesday, November 14, 1990 
1.00 PM 

Building 5, Room 119 

Calculations of complex shear layers give the best results when the con- 
stant in Rotta’s return-to-isotropy model is given values greater than those 
deduced from the measurements in grid turbulence. New data from grid 
turbulence are presented and assessed in the light of previous experiments. 
The results show that homogeneous turbulence decay is associated with 
increasing Reynolds number. It is argued that this finding has important 
implications for Reynolds-stress Modelling. Recent theoratical and experi- 
mental studies of the analogous model for pressure scrambling in the scalar- 
flux equations reveal even larger discrepancies in them “Monin constant” . 
The implications for the future of second-moment modelling are discussed. 


For additional information contact: T.-H. Shih, PABX 3-6680 


41 



INSTITUTE FOR COMPUTATIONAL 
MECHANICS IN PROPULSION 


ICOMP SEMINAR SERIES 


CURRENT TRENDS IN TURBULENCE MODELLING 

by 

W. Rodi 

University of Karlsruhe 
Institute for Hydromechanics 

Thursday, November 29, 1990 
1 : 00 - 2:00 p.m. 

ERB Building 5, Room 119 


A brief review is given of recent work in the area of modelling turbulence in near-wall regions and by 
Reynolds-stress-equation models. Various low-Reynolds-number versions of the k-e model and their 
damping functions are examined with the aid of results from direct numerical simulations and they are 
compared with respect to their performance in calculating boundary layers under adverse and favorable 
pressure gradients. A two-layer model is presented in which near-wall regions are resolved with a one- 
equation model and the core region with the standard k-e model. Various applications of this model 
are shown. The ability of the various models to simulate laminar-turbulent transition in boundary layers 
is discussed. Recent applications of a basic Reynolds-stress-equation model to complex flows of 
practical interest are presented. Finally, some recent proposals for improved Reynolds-stress-equation 
models are outlined and an outlook on possible future turbulence model developments is given. 


Contact: Charles Feiler, PABX 3-6681 
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COMPUTATION OF TRANSITIONAL AND TURBULENT 
FLOWS IN COMPLEX GEOMETRIES 


George Em Karniadakis 


Mechanical and Aerospace Engineering, 
Program in Applied and Computational Mathematics, 
Princeton University, Princeton, NJ. 


Friday, March 22, 1991 
2:00 - 3:00 PM 
Building 5, Room 119 


General, high-order numerical schemes are formulated appropriate for simulation of tran- 
sitional and turbulent incompressible flows in complex geometries. In particular, the new 
schemes are based on mixed explicit/implicit stiffly stable time-stepping methods and ex- 
plicit treatment of the pressure boundary condition that lead to enhanced stability and 
arbitrarily high-order time- accuracy dictated entirely by the employed integration rule. Hy- 
brid spectral element methods are then used for the spatial discretization of the variable 
properties governing equations in three-space dimensions. Special Neumann/ viscous sponge 
type boundary conditions are developed for open flows. Large or subgrid scales in the high 
Reynolds number regime are modeled through renormalization group (RNG) techniques. 

Simulations are than performed to study the transitional and fully turbulent stages of 
spatially developing flows. Here, we consider the flow over a backward- facing step, and the 
three-dimensional flow past a circular cylinder. For the first flow, the secondary instability is 
first investigated and the three-dimensional transitional states are computed through direct 
simulation (DNS). Transport and subgrid RNG models are then employed to simulate the 
high Reynolds number flow and heat transfer. Comparisons with experimental data in both 
regimes are presented. For the second, three-dimensional equilibria are accurately resolved 
via DNS; the series of bifurcations followed is simulated until the cylinder wake becomes 
turbulent. Our results suggest a successive period doubling in the temporal response of the 
flow, which eventually leads to a disordered state and renders the wake turbulent. 


For additional information contact: T.-H. Shih, PABX 3-66S0 
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A DYNAMIC MODEL FOR LARGE EDDY SIMULATION 
OF TURBULENT FLOWS 

Parviz Moin 

Stanford University 
and NASA Ames Research Center 



Tuesday, March 26, 1991 
10:00 - 11:00 AM 
Room 176, Sverdrup Building 


One major drawback of the subgrid scale models currently used in large 
eddy simulations is their inability to correctly represent with a single uni- 
versal constant different turbulent flows in rotating or sheared flows, near 
solid walls or in transitional regimes. A new subgrid scale model will be 
presented which alleviates many of these drawbacks. The model coefficient 
is inputted dynamically rather than inputted a priori. The basic idea in 
the derivation of the model is the utilization of the spectral information 
which is computed directly. This rather rich spectral information is not 
available in methods based on Reynolds averaged equations. The subgrid 
scale stresses obtained using the proposed model vanish in laminar flow and 
at a solid boundary. The results of large eddy simulation of transitional and 
turbulent channel flow that use the proposed model are in good agreement 
with the direct numerical simulation data. The same model was applied 
to the decay of isotropic turbulence with excellent agreement with the ex- 
perimental data. The model has been extended to compressible flows. A 
dynamic subgrid scale turbulent Prandtl number was derived. Its depen- 
dence on molecular Prandtl number, direction of scalar gradient, and the 
distance from the wall are in accordance with direct simulation. 


For additional information contact: T.-H. Shih, PABX 3-6680 
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TURBULENCE AND DETERMINISTIC CHAOS 

Robert G. Deissler 

Lewis Research Academy 
NASA Lewis Research Center 


Wednesday, April 10, 1991 
3:00 - 4:00 PM 
Room 119, Building 5 


Several turbulent and nonturbulent solutions of the Navier-Stokes equations are obtained. 
The unaveraged equations are used numerically in conjunction with tools and concepts from 
nonlinear dynamics, including time series, phase portraits, Poincare sections, largest Lia- 
punov exponents, power spectra, and strange attractors. 

Initially neighboring solutions for a low-Reynolds-number fully developed turbulence are 
compared. The turbulence, which is fully resolved, is sustained by a nonrandom time- 
independent external force. The solutions, on the average, separate exponentially with time, 
having a positive Liapunov exponent. Thus, the turbulence is characterized as chaotic. 

In a search for solutions which contrast with the turbulent ones, the Reynolds number (or 
strength of the forcing) is reduced. Several qualitatively different flows are noted. These are, 
respectively, fully chaotic, complex periodic, weakly chaotic, simple periodic, and fixed-point. 
Of these we classify only the fully chaotic flows as turbulent. Those flows have both a positive 
Liapunov exponent and Poincare sections without pattern. By contrast, the weakly chaotic 
flows, although having positive Liapunov exponents, have some pattern in their Poincare 
sections. The fixed-point and periodic flows are nonturbulent, since turbulence, as generally 
understood, is both time-dependent and aperiodic. 


For additional information contact: T.-II. Shih, PABX 3-6680 
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A CRITICAL ASSESSMENT OF TURBULENCE MODELING 


Charles G. Speziale 

ICASE, NASA Langley Research Center 


Monday, April 22, 1991 
1:30 - 2:30 PM 

Room 175, Sverdrup Building 


A variety of turbulence models including zero, one and two-equation mod- 
els as well as second-order closures will be reviewed. Based on comparisons 
with results from physical and numerical experiments on homogeneous tur- 
bulence, a strong case will be made for the superior predictive capabilities 
of second-order closure models. It will be shown how some significant im- 
provements in second order closure models have been recently achieved 
by means of invariance arguments coupled with a dynamical systems ap- 
proach. Several applications will be discussed including recent extensions 
to high speed compressible flows that were developed in connection with 
the National Aerospace Plane Project. 

For additional information contact: T.-H. Shih, PABX 3-6680 
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A CRITICAL-LAYER THEORY 
FOR BOUNDARY-LAYER TRANSITION 

Reda R. Mankbadi 


Lewis Research Academy 
NASA Lewis Research Center 


Tuesday, April 30, 1991 
10:00 - 11:00 AM 
Room 119, Building 5 

An asymptotic critical- layer theory is developed to study nonlinear interactions of a triad 
of instability waves leading to boundary- layer transition. This triad consists of a spatially 
growing plane fundamental wave and a pair of symmetrical, subharmonic oblique waves. The 
spatial development of the waves is determined by nonlinear viscous flow in the critical layer. 
The theory successfully captures not only the linear and parametric resonance stages, but also 
the later fully interactive regime of the transition process including the saturation and decay 
stages. The analysis is fully nonlinear, in that all the nonlinearly generated waves that were 
not originally present, are accounted for in the analysis. The three-dimensional nonlinear 
modifications of the mean flow are also considered. The theory applies to both ribbon- 
induced and naturally occurring transition. The analytically obtained amplitude equations 
are highly accurate but simple enough to be used for practical transition predictions. Results 
presented explain experimental observations and reveal novel features of the phenomena. 


For additional information contact: T.-H. Shih, PABX 3-6680 
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Research Briefs - 1990 


B.2 CMOTT Technical Meetings 

These axe CMOTT’s biweekly group meetings. In each meeting, one 
speaker from either CMOTT or other local modeling efforts at the NASA 
LeRc presents new ideas or work under progress. The informal seminars were 
intended not noly to keep the members informed of the latest development 
of local turbulence and transition modeling research but also to increase 
interactions between group members and other researchers at the NASA 
LeRc. 

The schedule of these biweekly events were widely disseminated before 
each series began. The following is the meeting schedule during the reporting 
period. 


48 


CENTER FOR MODELING OF 
TURBULENCE AND TRANSITION 

December 19, 1990 

CMOTT Members and SVR and IFMD Staff 
William W. Lion (6682) 

CMOTT Biweekly Meeting 

The following is a tentative schedule for CMOTT biweekly get-together from Dec. 
19, 1990 to Feb. 27, 1991. In each meeting, one member of our group will give an informal 
presentation on the subject she/he is interested in or currently working on. This will not 
only give us opportunities to interact with each other, but als'o keep us informed about 
exciting approaches in modeling turbulence other than our own specialities. However, the 
format of the meeting after the current period is open to suggestions. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


Dec. 19, 1990 

Zhigang Yang (61-2123) 

Understanding RNG in Turbulence: A Preliminary 
Report 

Jan. 2, 1991 

Tsan-Hsing Shih (6680) 

Realizability Concept and Its Applications in 
Turbulence Modeling 

Jan. 16, 1991 

Wai-Ming To (5937) 
Spectral Methods 

Jan. 30, 1991 

Le Tran (61-6701) 

Stagnation Point Turbulence and Heat Transfer 

Feb. 13, 1991 

William W. Liou (6682) 

Weakly Nonlinear Models for Turbulent Free Shear Flows 

Feb. 27, 1991 

Chris Steffen (8508) 

Optimizing Laminar Accuracy 
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Date: 

To: 

From: 


CENTER FOR MODELING OF 
TURBULENCE AND TRANSITION 

March 13, 1991 

CMOTT Members and SVR and IFMD Staff 
William W. Liou (6682) 


Subject: CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from 
March 20, 1991 to May 1, 1991. 

This round of biweekly meetings will be more of presentations of progress and hurdles 
in pursuing the subjects than presentations of results. Therefore, active participation is 
expected from the attendants. Note that these meetings are different from the CMOTT 
Seminar Series, whcih are mainly formal presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
are welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will staurt at 4:00 p.m. in Room 228, Sverdrup Building. 


Mar. 20, 1991 K. Zaunan (5888) 

Low Frequency Fluctuations in Separated Flows 
- Experimental Evidence 


Apr. 3, 1991 K. Kirtley (61-6659) 

Let’s Talk RNG 

Apr. 17, 1991 T.-H. Shih (6680) 

Lumley’s New Formulation of Dissipation Equation 


May 1, 1991 K. Kao (5965) 

Stability of Compressible Couette Flows 
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Date: 

April 25, 1991 

To: 

CMOTT Members and SVR and IFMD Staff 

From: 

William W. Liou (6682) 

Subject: 

CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from 
May 1, 1991 to June 12, 1991. Please note that Dr. Sliih’s and Dr. Kao’s talks have been 
rescheduled. 

The presentations will be informal and active participation is expected from the atten- 
dants. These meetings complement the CMOTT Seminar Series, which are mainly formal 
presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
are welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


May 1, 1991 T.-H. Shih (6680) 

Lurnley’s New Formulation of Dissipation Equation 

May 15, 1991 A. Shabbir (5927) 

On Some Closure Assumptions Regarding Time Scales, 
Local Equilibrium and Pressure Diffusion 


May 29, 1991 Z. Yang (61-2123) 

Near Wall k — e Modeling: Another Crusade 

June 12, 1991 K. Kao (5965) 

Stability of Compressible Couette Flows 
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THE STUDY OF PDF TURBULENCE MODELS 

IN COMBUSTION 


Andrew T. Hsu 

Sverdrup Technology, Inc. 

& Center for Turbulence Modeling 
NASA Lewis Research Center 
Cleveland, Ohio 


9th National Aero-Space Plane 
Technology Symposium 
November 1-2, 1990 


Paper Number 107 
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I. Introduction 


(U) In combustion computations, it is well known that the predictions 
of chemical reaction rates (the source terms in the species conservation equa- 
tions) are poor if conventional turbulence models are used. The main dif- 
ficulty lies in the fact that the reaction rate is highly non-linear, and the 
use of averaged temperature, pressure and density produces excessively large 
errors. Moment closure models for the source terms have attained only lim- 
ited success. The probability density function (pdf) method seems to be 
the only alternative at the present time that uses local instantaneous val-. 
ues of the temperature, density, etc., in predicting chemical reaction rates, 
and thus is the only viable approach for more accurate turbulent combustion 
calculations. 

(U) The fact that the pdf equation has a very large dimensionality ren- 
ders finite difference schemes extremely demanding on computer memories 
and thus impractical, if not entirely impossible. A logical alternative is the 
Monte Carlo scheme, which has been used extensively in statistical physics. 
The evolution equations for the joint pdf of the velocity and species mass 
fraction have been successfully solved using Monte Carlo schemes, see, e.g., 
Pope[l]. However, since CFD has reached a certain degree of maturity as 
well as acceptance, it seems, at least from the stand point of practical appli- 
cations. that the use of a combined CFD and Monte Carlo scheme is more 
beneficial. Therefore, in the present study a scheme is chosen that uses a 
conventional CFD flow solver in calculating the flowfield properties such as 
velocity, pressure, etc., while the chemical reaction part is solved using a 
Monte Carlo scheme. 

(U) A combined CFD- Monte Carlo computer algorithm has been de- 
veloped recently. As a first calibration of the Monte Carlo solver recently 
developed in this work, the discharge of a heated turbulent plane jet into 
quiescent air is studied. Experimental data for this problem shows that when 
the temperature difference between the jet and the surrounding air is small, 
buoyancy effects can be neglected and the temperature can be treated as a 
passive scalar. The fact that jet flows have a self-similar solution lends con- 
venience in the modeling study. Furthermore, the existence of experimental 
data for turbulent shear stress and temperature variance (temperature fluc- 
tuation) make the case ideal for the testing of pdf models wherein these 
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values can be directly evaluated. 

(U) The following section presents the methodologies used in this study, 
and a discussion on the numerical results and their comparison with experi- 
mental data follows. 


II. Method 

2.1 Governing equations 

(U) The Favre averaged momentum and species transport equations can 
be written as 

pd t Ui + pUjdjUi = -dip + vdjTij — dju'-u" 

pd t Yi + pujdjYi = pdj{DdjY t - Y"u") + dj(DdjpdjY") + pw x 

where u t is the velocity and Y x is the mass fraction. The major problem en- 
countered in the numerical study of turbulent combustion lies in the modeling 
of the chemical source term. pw ix known as the chemistry closure problem. 
The Source term W{ is an exponential function of temperature, and it is well 
know that the use of averaged temperature, T, in evaluating pW{ can cause 
egregious errors. Therefore the accurate prediction of turbulent combustion 
using conventional turbulence models becomes very difficult. This motivated 
the use of pdf methods. If the pdf, P. is given, then the mean of the source 
term can be evaluated exactly: 

pw, = J...Jpw t (Y\ Y n ,T,p)P(Y l ,...,Y n ,T,p)dY 1 ...dY n dTdp 

To solve for the pdf, we need the following: 

2.2 Evolution equation for the pdf 

(U) The evolution equation for the probability density function of the 
mass fractions, temperature, and density, P(Y\, />), can be written 

as[l] 

;V 

pd t P + pv a d Q P + pJ2^. Pn)P) 

l = l 

= -d*(p < > p) - p > P) 

*=i j=i 
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where the terms represent mean convection, chemical reactions, turbulent 
convection, and molecular mixing, respectively; P is the density-weighted 
joint pdf: 

P = pP , 

6 is the scalar dissipation: 

t{j pD{jdj\'{ , 

0 t -’s are variables such as the mass fractions, density, and temperature, and 
<, > denotes the mathematical expectation of a function. 

(U) The left hand side of the above equation can be evaluated exactly 
and requires no modeling; the right hand side terms contain the conditional 
average of the Favre velocity fluctuation and the conditional average of the 
scalar dissipation and require modeling. 

2.3 Finite difference solution for the NS equations 

(U) In order to simplify the problem so a s to concentrate on the study of 
the pdf models, we have confined our numerical procedures to parabolic flows. 
A parabolic NS solver with a k — t turbulence model has been developed. 
The k — t model has been tested by solving standard cases such as flat plate 
boundary layers and free shear layers with satisfactory numerical results. 

2.4 Monte Carlo scheme for the pdf equation 

(U) A grid dependent Monte Carlo scheme has been employed, primarily 
following Chen and Kollmann[5]. The task here is to construct an ensemble 
of sample points, each sample has its own distinct properties such as tem- 
perature. mass fraction, etc.: these properties change with time or location 
such that the probability function of the ensemble evolves according to the 
evolution equation. Consequently, the pdf of the ensemble is the desired 
approximate solution for the pdf equation. 

(U) Suppose N samples are assigned to a grid cell in the flow domain, 
the pdf for the ensemble of N samples can then be written as 

= 7 ? 

1 j= i 

where o 3 is the scalar function value carried by the j th sample, e.g., the mass 
fractions, etc. 
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(U) In order to have the pdf P m evolve according to the pdf evolution 
equation, we discretize the equation on a given grid and write, for parabolic 
flows: 

P m x+dx,j = &jP*x.j + l + 3jP* X 'j + 'Vj P m xj — l 

and require 

+ Pj + a = i 


2.5 Recontamination in a Monte Carlo scheme 

(U) Using a very simple test case of a convection/diffusion process with 
two scalars, it was found that the previous scheme[5] does not conserve mass 
fractions due to re-contamination. It is found that in order to conserve the 
mass fractions absolutely, one needs to add further restriction to the scheme, 
namely 

a, + 7j = <*j- 1 + 7j+i (*) 

A new computer algorithm was devised and tested. This algorithm uses a 
few extra arrays to store informations form the previous time level and thus 
eliminates the repetition in the sampling process. With the simple test case 
of two scalars with assumed constant coefficients in the pdf equation, the 
new algorithm is shown to conserve the mass fractions perfectly. Deficiencies 
such as directional bias and re-contamination that were found in the previous 
algorithm are completely eliminated. 

(U) It is not yet known whether one can indeed devise a scheme that 
satisfies relation (*) in general flows, where the coefficients of the pdf equation 
are variable and are calculated from the flow velocities, turbulence time scale, 
etc. 
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III. Results 


(U) The sketch of a heated plane jet is given on this page. The flow 
domain is divided into 30 cells in the cross direction, and 100 samples are 
assigned to each cell. The numerical results from the present study are 
compared with experimental data in the following figures. 



Figure 1. Sketch of a heated turbulent free jet. 
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(U) The numerical solution for the ensemble averaged pdf’s of the tem- 
perature: y/yo.s — 0 - where y 0 . 5 is the location where T — 0.5To . In order to 
interpret the data, the ensemble averaged pdfs are integrated to obtain the 

approximate continuous pdf: P = P m (AT)/ AT. The result of this integra- 
tion is shown as vertical bars. The normal distribution is plotted as a dotted 
line for comparison. 




Fig. 2 probability density distribution 
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(U) The numerical solution for the ensemble averaged pdf’s of the tem- 
perature: y/yo.s = 0.44, where y 0 . 5 is the location where T = 0.5T o . In order 
to interpret the data, the ensemble averaged pdfs are integrated to obtain 
the approximate continuous pdf: P = P m (&T)/ AT. The result of this in- 
tegration is shown as vertical bars. The normal distribution is plotted as a 
dotted line for comparison. 


Fig. 3 probability density distribution 



(T-<T>)/a 
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(U) The computed velocity profile is compared with experimental data 
from Foerthmannf?]. The good agreement between the numerical solution 
and experiments indicates that the parabolic flow solver employed in the 
present study is fairly accurate and can serve its purpose. 


• Fig. 4 Velocity distribution 
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(U) As a calibration tor the k — e turbulence model, the turbulent 
shear stress from the present study is compared with measurements given 
by Brandbury[S], and by Gutmark and Wvgnanski[9]. Except for a slight 
over prediction towards the edge of the jet, the numerical solution generally 
agrees well with the experimental data. This comparison is very crucial for 
pdf calculations because it is an indicator of the accuracy of the numerically 
predicted turbulent kinetic energy, k. and dissipation, e, which are used to 
determine turbulent time scale, t. from t = k/e. An accurate prediction for 
k and e is essential for the correct modeling of molecular mixing in the pdf 
calculation. 


Fig. 5 Turbulent shear stress 
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(U) The results for mean temperature from the Monte Carlo solution of 
the pdf equation for the temperature compared with experimental data from 
Ref. 10-12. 


'Fig. 6 Mean temperature 
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(U) The standard variance of the temperature, or the temperature fluctu- 
ation, compared with experimental data from Ref. 12 and 13. Because of the 
limited number of sample points (N=100) used in the computation, the nu- 
merical results are scattered: however, these results agree very well with the 
experimental observations. Particularly important is the agreement between 
the predicted and the measured temperature fluctuations, for in case of finite 
rate chemistry calculations, the instantaneous temperature will be used in 
determining the reaction rate, and the correct prediction of the temperature 
distribution will ensure the accuracy in reaction rate calculations. 


O 



y/y.5T 0 


64 


Acknowledgement 


The author wish to thank Prof. S.B. Pope of Cornell University and Dr. 
J.-Y. Chen of Sandia National Labs for their help. This work is supported 
by the NASA Lewis Research Center under Contract NAS3-25266. 

References 


1. Pope, S.P., Prog. Energy Combst. Sci. 1985. Vol.ll. pp. 119-192 

2. Kolhnann, W., Theoret. Comput. Fluid Dynamics, 90 1. pp. 249-285 

3. Chen, J.-Y. and Kollmann, W. Combust. Sci. and Tech. 1989, Vol. 
64, pp. 315-346 

4. Shuen, J.-S. and Yoon, S. A1AA J. Vol. 27 No. 12, 1989, pp. 1752-1760 

5. Chen, J.-Y. and Kollmann, W. Combust, and Flame. 1990, pp. 75-99 

6. Brabbs, T.A. and Musiak, J.D., Private Communications 

7. Schlichting, IT, “Boundary Layer Theory.” McGraw-Hill, 1989. 

8. Brandbury, L.J.S., JFM, Vol. 23, 1965, pp. 31-64. 

9. Gutmark, E. and Wygnanski, I., JFM, Vol. 73, 1976, pp. 465-495. 

10. Browne, L.J.S., Antonia, R.A., and Chambers, , A.J., JFM, Vol. 149, 
1984, pp. 355-373. 

11. Uberoi, M.S. and Singh, P.I., Physics of Fluids, Vol. 18, No. 17, 1975, 
pp 764-769. 

12. Jenkins, P.E. ASME J. of Eng. for Power., Vol. 98, 1976, pp.501-505. 

13. Antonia, R.A., et al., Int. J. of Heat Mass Transfer, Vol. 26, No. 1, 
1983, pp. 41-48. 


65 


On Recontamination and Directional-Bias Problems 
in Monte Carlo Simulation of PDF Turbulence Models 


Andrew T. Hsu 

NASA Lewis Research Center, Cleveland, Ohio 44135 
Phone 218/826-6648 


Turbulent combustion can not be simulated adequately by conventional moment closure turbulence models. 
The difficulty lies in the fact that the reaction rate is in general an exponential function of the temperature, and 
the higher order correlations in the conventional moment closure models of the chemical source term can not be 
neglected, making the applications of such models impractical. The probability density function (pdf) method 
offers an attractive alternative: in a pdf model, the chemical source terms are closed and do not require additional 
models. 

The partial differential equation for the probability density function, P, can be written as 

N 

pd t P + pv Q d Q P + P ^n)P) 

»=i 
N * 

. = -9a{p < Vo\A > P) > P) 

where the terms represent the time derivative, mean convection, chemical reaction, turbulent convection, and 
molecular mixing, respectively. The fact that the pdf equation has a very large dimensionality renders finite 
difference schemes extremely demanding on computer memory and CPU time and thus impractical, if not entirely 
impossible. A logical alternative is the Monte Carlo scheme, wherein the number of computer operations increases 
only linearly with the increase of number of independent variables, as compared to the exponential increase in a 
conventional finite difference scheme. 

A grid dependent Monte Carlo scheme following that of J.Y. Chen and W. Kollmann has been studied in the 
present work. In dealing with the convection and diffusion of the pdf, the pdf equation is discretized on a given 
grid, e.g., 

P*+dxJ = Q jP*j+l PjPxJ "b IjPzj-l 

where 

Qj + Pj + 7j = 1 

However, if this is the*only restriction satisfied by the numerical algorithm, the mass fractions may not be con- 
served due to re-contamination, and directional- bias also appears. These phenomena are illustrated in Figure 1: 
Consider a mixing layer; use white balls to represent contaminants in the upper stream and black balls to represent 
contaminants in the lower stream. As the two streams move toward right, the location of the white balls and black 
balls are interchanged randomly to simulate convection and diffusion. From Figure 1, it is clear that directional- 
bias caused by recontamination caused the center of the mixing layer to drift downward. ( DirectionaJ-bias can 
be partially corrected by changing sweeping directions. ) One also notices that after the first marching step, the 
conservation law is violated, reflected in the Figure as missing white or black balls. 

It is found that in order to conserve the mass fractions absolutely, one needs to add further restriction to the 
scheme, namely 

Q j + 7i = a ;-l + 7>+i 

A new algorithm was devised that satisfies this restriction in the case of pure diffusion or uniform flow problems. 
Using the same example, it is shown that absolute conservation can be achieved. This result is shown in Figure 2. 
One can see that the diffusion process is symmetric, and the problem of directional-bias is eliminated. 

Although for non-uniform flows absolute conservation seems impossible, the present scheme has reduced the 
error considerably. 
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Figure 1. Illustration of directional-bias and recontaui nation 
problems in Monte Carlo simulation. 

Y* 



Figure 2. Solution obtained by using additional constraint 
in the solution procedure. 
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PROGRESS IN THE DEVELOPMENT OF 
PDF TURBULENCE MODELS FOR COMBUSTION 


Andrew T. Hsu 

Sverdrup Technology, Inc. 
NASA Lewis Research Center 
Cleveland, Ohio 


10th National Aero-Space Plane 
Technology Symposium 
April 23-26, 1991 


Paper Number 213 


L Introduction 


(U) A combined Monte Carlo — CFD algorithm has been developed re- 
cently at NASA lewis Research Center for turbulent reacting flows. In this al- 
gorithm, conventional CFD schemes are employed to obtain the velocity field 
and other velocity related turbulent quantities, and a Monte Carlo scheme 
is used to solve the evolution equation for the probability density function 
(pdf) of species mass fraction and temperature. 

(U) In combustion computations, the predictions of chemical reaction 
rates (the source terms in the species conservation equations) are poor if 
conventional turbulence models are used. The main difficulty lies in the fact 
that the reaction rate is highly non-linear, and the use of averaged tempera- 
ture produces excessively large errors. Moment closure models for the source 
terms have attained only limited success. The probability density function 
(pdf) method seems to be the only alternative at the present time that uses 
local instantaneous values of the temperature, density, etc., in predicting 
chemical reaction rates, and thus may be the only viable approach for more 
accurate turbulent combustion calculations. The closure problem and the 
need for pdf has been discussed in detail in a previous paper[l]. 

(U) Assumed pdf are useful in simple problems; however, for more gen- 
eral combustion problems, the solution of a evolution equation for the pdf is 
necessary. 


IL Method 

(U) Conventional CFD flow solvers are quite common and needed no 
explanation here. We will concentrate on the derivation and solution of the 
pdf evolution equation in this section. 

2.1 PDF evolution equation of mass fraction. 

(U) For simplicity, the pdf equation for a single scalar will be derived 
here. The extension to multi species is trivial. The species transport equa- 
tions ( for a single species ) can be written as 

pd t Y = - pu k d k Y + pdk(DdkY) + pw (0.1) 

where u * is the velocity, Y is the mass fraction, and w is the chemical source 
term; summation for repeated indices is understood. 
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Let P(Y) be the probability density function, and define <f> = exp(z’Ac), 
then the characteristic function of the pdf can be written as the ensemble 
average of <f>: 

<t>=J P(Y)<f>dY , 

and the pdf can be written as the inverse Fourier transform of the character- 
istic function: 

P(Y ) = 2 “ / < ^ > exp{-i\Y)dY. 

Now differentiate <f> with respect to time, 

d t 4> = tAeip(t'Ar)(5 t r), 

and replace d t Y in the above equation with the right hand side of eq. (1), 
we obtain 

pd t <f> = i\exp{i\Y)[-pv. k d k Y + d k (pDd k Y ) + pw], 

which can be rearranged as 

dt(p<£) + d k (pu k <t>) = i\<f>d k (pDd k Y ) + i\<f>pw 

Take ensemble average of the above equation gives 

dt < p<f> > +d k < pu k (f> >= i\ < pwcf) > +t'A < <f>d k (pDd k Y ) > . 

The inverse Fourier transform of the above equation gives the evolution equa- 
tion for the probability density function: 

d t (pP) + d k (p < u k > P)+ dy(pwP) 

= -d k (p < u' k \Y > P)-M< d k ( P Dd k Y)\Y > P ) 

where the terms represent mean convection, chemical reactions, turbulent 
convection, and molecular mixing, respectively. 

For multi species and temperature, the pdf is written as P(ipi, ..., V\»)> 
where t/\', t = represent the scalars including mass fractions and tem- 

perature. A similar derivation gives 

N 

pdtP + pv a d a P + ...,xp N )P} 

•= 1 

= -d a (p < >p)~ > P )• 

.=1 j = 1 
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(U) The left hand side of the above equation can be evaluated exactly 
and requires no modeling; the right hand side terms contain the conditional 
average of the Favre velocity fluctuation and the conditional average of the 
scalar dissipation and require modeling. 

2.2 Modeled pdf equation 

(U) The first term that needs modeling is the turbulent diffusion term, 
for which the following gradient model is used: 


- < u' k \xj>i >P* D t d Q P , 

where D t is the turbulent diffusion coefficient, which is set to be equal to the 
eddy viscosity, i.e., the turbulent Schmidt number is equal to one. 

(U) The next term needs modeling is the molecular mixing term. A 
coalescence/ dispersion model is used for this term, which has the following 
general form: 

Cijilfe > p ) 

»= i i=i 

-^7 \ j p W)p{V)Tm\rwdr -p 

= D m P 

where T is the transition probability. For more detail description of the 
molecular mixing model, see Ref 2. The modeled the pdf equation is then 

N 

pdtP + pv a d a P + p ]T d^{wi(rp i, rp N )P} 

= d j (D t d a P) + D m P 


2.3 Solution of the pdf equation 

(U) A fractional step method is used to solve the pdf equation in a 
Monte Carlo simulation. Three processes, namely, (1) convection/diffusion, 
(2) molecular mixing, and (3) chemical reaction, axe simulated consecutively: 

N 

(dt + p'52d+ i {w i (tJ>i,...,iJ> N )})(d t -D m )(d t + pvd a - d a {D t d a )P = 0. 

1=1 
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This equation is solved using a Monte Caxlo simulation, i.e., the pdf is 
represented by an ensemble of samples. These sample points move in the 
hyperspace of (x 1 ,x 2 ,x 3 , V»i, ... ,V>„) according to the above equation. 


HI. Results 

(U) A heated turbulent plane jet had been calculated and results re- 
ported at last NASP symposium. The case has since been recalculated using 
modified turbulence model, and improved results are shown here. 




Turbulent shear stress 
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(U-U1 )/(U2-U t ) 


(U) A hydrogen-fluorine diffusion flame was calculated. The flowfield 
and flow conditions, as well as the calculated velocity and turbulent kinetic 
energy distribution, are given on this page. 
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(U) Figures show the calculated heat release of the H2-F2 flame (tem- 
perature distribution) compared to experimental data, and the species mass 
fraction from pdf calculation. 


Fig. 1 Temperature distribution: H2+F2=2HF 
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(U) Computed pdf of temperature distribution in the H2-F2 flame. 




a^/a T ad , 
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A STUDY OF HYDROGEN DIFFUSION FLAMES USING 
PDF TURBULENCE MODEL 


Andrew T. Hsu* 

Sverdrup Technology, Inc. 

NASA Lewis Research Center, Cleveland, Ohio 44145 


Abstract 

The application of probability density function (pdf) turbulence models is addressed in this work. For 
the purpose of accurate prediction of turbulent combustion, an algorithm that combines a conventional CFD 
flow solver with the Monte Carlo simulation of the pdf evolution equation has been developed. The algorithm 
has been validated using experimental data for a heated turbulent plane jet. The study of H 2 -F 2 diffusion 
flames has been carried out using this algorithm. Numerical results compared favorably with experimental 
data. The computations show that the flame center shifts as the equivalence ratio changes, and that for the 
same equivalence ratio, similarity solutions for flames exist. 


Introduction 


It is the consensus of the combustion profession 
that the prediction of chemical reaction rate (the 
source term in a species conservation equation) is 
poor if a conventional turbulence model is used. The 
main difficulty lies in the fact that the reaction rate 
is highly non-linear, and the use of averaged tem- 
perature, pressure and density produces excessively 
large errors. Moment closure models for the source 
terms have attained only. limited success because the 
assumptions for such models to be valid often can 
not be satisfied. The probability density function 
(pdf) method seems to be the only alternative at the 
present time that uses local instantaneous values of 
the temperature, density, etc., in predicting the chem- 
ical reaction rate, and thus is the only viable approach 
for more accurate turbulent combustion calculations. 
Two main lines are being followed in pdf methods: 
one uses an assumed shape for the pdf, the second 
solves a pdf evolution equation; the present paper ad- 
dresses only the latter. There has been significant 
progress in the study of pdf turbulence models in low 
speed flows in the past decades. These developments 
were summarized in Refs. [1], [2] and [3]. In spite 
of this progress, pdf turbulence modeling remains a 
nascent discipline with many unresolved issues. 

The fact that the pdf equation has a very large 
dimensionality renders the use of finite difference 
schemes extremely demanding on computer memories 
and thus impractical, if not entirely impossible. A log- 
ical alternative is the Monte Carlo scheme, which has 
been used extensively in statistical physics. The evo- 
lution equations for the joint pdf of the velocity and 
species mass fraction have been successfully solved 
using Monte Carlo schemes, see, e.g., Pope[l]. How- 
ever, since CFD has reached a certain degree of ma- 
turity as well as acceptance, it seems, at least from 
the standpoint of practical applications, that the use 
of a combined CFD and Monte Carlo scheme is more 
beneficial. Therefore, in the present study a scheme 
is chosen that uses a conventional CFD algorithm to 
solve the Navier-Stokes equations and provide flow- 
field properties such as velocity, pressure, etc., and the 
chemical reactions are calculated by using a Monte 


Carlo scheme to solve a pdf evolution equation. 

The combined CFD-pdf solver has been devel- 
oped recently and validated using non-reacting flow 
data for a heated turbulent plane jet. The algorithm 
has been further tested in the numerical study of an 
H 2 -F 2 diffusion flame. This diffusion flame was stud- 
ied experimentally by Mungal &; Dimotakis[4] and 
Hermanson Dimotakis[5]. The numerical results 
from the present study are compared with these ex- 
perimental data. 


Theory 

Governing equations for reacting flows. 

Flows with chemical reaction are governed by 
the continuity equation, momentum equations, en- 
ergy equation, and species transport equations: 


dtp 4- dipui = 0 

pd t U{ -j- pujdjUi — -d x p 4- pdjTij 
pd t h 4- puj dj - d t p- Uj dj p 

— —djqj -F Q 

pd t Yk 4- pUjdjYk = pdj(DdjYk) + w k 

k = 1 , 2 , 


where u, is the velocity, Y k is the mass fraction, and 
Wk is the chemical source term. (In addition to these 
equations, one also needs the equation of state.) For 
turbulent flows, we substitute 

U { = Ui 4- u • , 

V. = Yi + vy, 

into the above equation and take an ensemble aver- 
age. In the process of averaging, new unknown quan- 
tities in the form of correlations appear, e.g., 
ujyy, etc. These quantities can all be modeled using 
conventional turbulence models, such as two equation 
models or second order closure models. A problem 
unique to flows with chemical reaction is the average 
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of the chemical source term, pWk , The major diffi- 
culty lies in the fact that Wk is a highly nonlinear 
(usually exponential) function of the temperature. It 
is well know that the use of averaged temperature, T, 
in evaluating puJi can cause egregious errors. One may 
consider the effect of the temperature fluctuation, T' , 
by applying a moment closure model to the term puJf; 
however, such a closure model results in an infinite se- 
ries that converges only when T a ~ T and T* *C T. In 
many combustion problems, these two conditions are 
violated: in fact, we often have T a T and T* ~ T 
instead. In view of the above, the prospect of an ac- 
curate prediction of turbulent combustion using con- 
ventional turbulence models seems dismal. The above 
fact motivated the use of pdf methods. 

Pdf evolution equation 

Given a set of m random vari- 
ables, •••,0m * and the joint probability den- 

sity function, P(rp\ t rp 2 , * * • , x, y, z, t), the mean of 
any random function, /(^i, V^i***, V’m;*, y, can 
be calculated as 

f(x,y,z,t) = J-'-J f Pdifiidfa ■■■ dxl> m . (3) 

For simple flows, one could assume the shape of 
the pdf and compute the source term, puT, based on 
the assumed pdf using the above integral. For more 
general problems, one needs to solve a pdf evolution 
equation for a more accurate pdf distribution. In the 
latter case, the evaluation of the source term, ~puJk, is 
no longer necessary because the mean values of the 
temperature, species mass fractions, which we previ- 
ously would obtain by solving the transport equations 
(Eq. 1), can now be evaluated directly from the pdf 
using the above integral. 

The pdf evolution equation can be derived from 
the transport equations (Eq. 1) in many different 
ways, using a Dirac delta function, a characteristic 
function, or a characteristic functional. [1-3] The evo- 
lution equation of a single point probability density 
function of scalar random variables rp 1 , • • • , \p m can be 
written as 

m 

pd t P + pv a d a P + p "^2 • •• , xp m )P) 

i=l 

= -d a (p < = *k>P) 

m m 

= d>t>p) ( 4 ) 

•=i j=i 

where the terms represent the rate of time change, 
mean convection, chemical reaction, turbulent con- 
vection. and molecular mixing, respectively; P is the 
density-weighted joint pdf: 

P = PP/P, 

tij is the scalar dissipation: 

Uj = Dd a <j>id a 4>j , 


(where D is the diffusion coefficient), and < x|y > 
denotes the mathematical expectation of a random 
function x conditioned upon y. 

The left hand side of eq. (4) can be evaluated 
exactly and requires no modeling; the right hand side 
terms contain the conditional expectation of the ve- 
locity fluctuation and the conditional expectation of 
the scalar dissipation, which are new unknowns and 
require modeling. 

Closure models for pdf equation 

The first term that needs modeling is the tur- 
bulent convection term, for which one can use the 
following gradient model: 

-<v n a \ rp k > P^DtdaP, 

where D t is the turbulent diffusion coefficient, which 
is set to be equal to the eddy viscosity, i.e., the tur- 
bulent Schmidt number is equal to one. 

The next term needs modeling is the molecular 
mixing term. A coalescence/dispersion model is used 
for this term, which has the following general form[6]: 

N N 
»=i j = i 

J J Pwpwmtf, vwdr - p ’ 

= M(P) 

where T is the transition probability. A new mix- 
ing model continuous in time is recently developed by 
Hsu and Chen[7]. For more detail description of the 
molecular mixing models, see Refs. 6 and 7. 

The modeled pdf equation is then 

m 

pd t P + pv a d a P + P^2d^ t {wi(rpi, ■ ■ ■ ,xp m )P} 

i=i 

=d j (D,d a P) + M(P), (5) 

which can be solved using a Monte Carlo simulation. 

Numerical Methods 

In the present study, a combined finite 
difference-Monte Carlo solver is developed. For the 
velocity field, the Navier-Stokes equations and a k — c 
turbulence model are solved using a finite difference 
method; for scalar variables such as the temperature, 
mass fractions, etc., the pdf evolution equation is 
solved using a Monte Carlo scheme. The CFD flow 
solver provides the Monte Carlo solver with the mean 
velocity and a turbulence time scale r, where r = Jb/e, 
and the Monte Carlo solver provides the mean flow 
solver with the density, p. 

The solution’of the Navier-Stokes equations 


78 


Since the primary concern of the present work is 
pdf modeling, we choose the simplest possible solver 
for the N-S equations. The equation is transformed 
into a general coordinate system using the following 
general relation: 

y = y(£, n) 

and 



The transformed momentum equation for steady 
flows is 

u d(U + (v - iiy^—drjU = — \d&+ —d^i/rd^u), 

yq P yq 

where the viscous term involves gradient in the £- 
direction is neglected since we are only interested in 
shear flows in the present study. 

First order upwind difference is used in the 
^-direction so that a marching scheme for steady 
parabolic flows can be used. To ensure stability, a 
flux splitting scheme is used in the 77 -direction. 

The solution of the pdf ecmation 

A fractional step Monte Carlo method, as will 
be described in the following, is used to solve the pdf 
equation. We first discretize the time derivative in eq. 
(5) using finite difference: d t P = (P n+I — P n )/At. 
then the pdf evolution equation (eq. 5) can be written 
as 

m 

P n+1 = {1 - A*D 0 3 0 - A l 

1=1 

+Atd 0 D,d a + AtM } P n 


(3) chemical reaction: 

P n+1 = (1 + AtR)P mm 

In the present study, only steady flows are considered, 
so d t is replaced by d( and a marching scheme in the 
^-direction is employed. 

In a Monte Carlo simulation, the continuous pdf 
is replaced by N x M delta functions, 

2 , • • • . V m ;x,y,z,t) 

n= 1 

X <5(V>2 - 4 n) (0) ■ • • Wm - 

where each product of the m delta functions repre- 
sents one event of an ensemble of N sample events. 
An event can be thought of as a fluid particle, and 
the evolution of P* entails the movement of the par- 
ticles in the physical space as well as the phase space 
(t/> — space). The movement of the particles is, of 
course, governed by the pdf evolution equation, and 
is simulated in the following three steps. 

Step 1: convection 

Replace P by P m in the pdf evolution equation 
and transform it into the £ — rj coordinate system, the 
convection process for a steady flow can be written as 

pud{P~ = -p(v - uy{)—dr)P m 4- —d tl {pD t dr 1 P m ) i 
yq yq 

Discretizing the above equation using a finite differ- 
ence scheme, we can write 


Using approximate factorization, the above equation 
is recast as 


Pi.j ~ Q Pi-l.j + l 


■3 PL 


1 ■} 


rPL 


i.>- 1 


P n+l =(l - At y a.,, wi(xl> i,---.0 m )) 

i=l 

x(l + Aa/)(1 - A tpvda T A td a D,d a ) P n 

+0( At) 

= (1 + A*C)(1 + AtM)(l + A tR)P n +0(At), 

where C denotes the convection operator, m the 
molecular mixing, and R chemical reactions. With 
the above expression, we can three processes consec- 
utively: 

(1) convection: 

P m = (1 + AfC)P n , 

(2) molecular mixing: 


(Here again we used a one sided difference so that a 
marching scheme could be used.) The above equation 
states that if we divide the flow field into cells, then 
the pdf at point (i,j) can be written as a linear com- 
bination of the pdf at neighboring points, To simulate 
this process with a Monte Carlo scheme, we move the 
sample particles between cells according to the above 
equation. For instance, the particles of cell (x, j) will 
be obtained by choosing randomly aiV particles form 
cell (x — 1 , j -f 1 ), pN particles from cell (x — 1 ,;), and 
7 N particles from cell (x — 1 ,; — 1 ). In order for the 
total number of particles not to change, we require 
that a + /? + 7 = 1 . 

It is worth noting that this method can easily 
be extended to elliptic flows and applied to general 
curvilinear coordinate systems. 


P mm = (1 + A tM)P\ 


Step 2: molecular mixing 
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The molecular mixing process is simulated by the 

following binary interaction model: 

pud{P 

P(^P , )P(r)T(lP\cP\r)d^P , d^l; , ' - P 

where T is the transition probability. By assigning 
various functions to T, we would have different mix- 
ing models. In the present study, the modified Curl 
model by Janicka et al. [8] is used. In this model, the 
transition probability is given as 


will acquire new composition at the next moment, and 
the changing rate is 

d<f>i 

= 1,2, • • • , M, 

where is the reaction rate for that specific sample 
point. If we regard u/ t - as the convection velocity of 
the particle in the phase space, the it is clear that 
the above ODE represent the movement of a sample 
particle in the phase space. 

The various chemistry models considered in this 
study are presented in the following section. 



T{* 


i 

0 


for xl>' < x/j < i/>'\ 
or t/>” < ip < 
otherwise. 


with Cd — 6.0. 

In the Monte Carlo simulation, the above model 
is realized in the following manner: Divide the flow 
domain into small cells, each containing N sample 
particles. Given a small time interval At and a tur- 
bulent time scale r, select randomly N mT pairs of par- 
ticles, where 


= O.S^Lv. 

Ct 


and let a pair, say, m and n. mix as follows 


6 n (t 4- At) = A<j> m (t) + (I - A)4> n (t) 
+ At) = A<f> n (t) + (1 - A)4>m{t) 


where .4 = 0.5£, with £ a random variable uniformly 
distributed on the interval [0.1]. The remaining N — 
2 Nmx particles remain unchanged: 

<f> n (t + AO = O „(0 


The turbulent time scale is supplied by the finite 
difference solution of a k — c turbulence model: r = 
k/e. 

This model admits the non-physical jump con- 
dition and does not produce the correct long time 
behavior for decay problems in homogeneous turbu- 
lence. A continuous model that predicts the correct 
long time behavior for turbulence decay problem has 
been introduced by Hsu and Chen [7]. But as shown 
in Ref. [7], for practical combustion problems where 
the long time statistical behavior is not crucial, the 
modified Curl model gives acceptable results. In the 
present study, both the modified Curl model and the 
continuous model have been used. 


Step 3: chemical reaction 

Chemical reaction is represented in the Monte 
Carlo simulation by the movement of sample particles 
in the phase space due to reaction. A sample point 
with a given composition {oi,<^ 2 i * * at time t 


Chemistry Models for //? — Fn Reaction 
Finite rate chemistry . 

The reaction of Hn + F 2 = 2 H F can be repre- 
sented by the following two step chain reactions[4]: 

H 2 + F — • kl HF+H 

H + F 2 — ^ HF + F 


with 


*1 = 2.6 x 10 12 r°' s exp , 


*2 = 3 X lO 9 ^ 5 exp 


— 1680 \ 
RT ) 


where T is in K, R in cal mol~ l K ~ l , and k in cm 3 
mol~ l sec -1 . 

The above expressions only provide the forward 
reaction rate; the backward reaction rate can be cal- 
culated using the equilibrium coefficient: 


kb kf / k C q , 


and the equilibrium coefficient can be evaluated using 
Gibbs free energy. Polynomials for Gibbs free energy 
for various species can be found in Ref. [8]. A nu- 
merical experiment shows that the backward reaction 
rates are much smaller than that of the forward re- 
action for the two reaction steps described above and 
can be neglected. 

Let C denote the mole concentration; for each 
sample particle, we need to solve the following set of 
ODEs: 

dCn 7 u n r 

- -k 1 Ch? Cf = Ui 


dt 

dC F , 

dt 
d Ch 
dt 
d C F 


— —k o Ch C F j — ujn 


— — -f- U) 2 — CJ 3 


dt 
dCurF. 

dt 


— <J| — CO 2 — CJ4 


= —to l — CO 2 — UJ 5 


dk o 

■ 3 T = -£ 

1=1 


where h is the. enthalpy and h ^ is the heat of forma- 
tion. 
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Using n to denote the time level, we found that 
the following semi-implicit scheme seems to be most 
stable: 


a few hundred degrees Kelvin. Under these circum- 
stances, the variation of specific heat C r is negligible, 
and the temperature can be calculated as 



C£' 


C]t 


1 


c? +1 


/^ri + 1 


r*2 

1 + AtkyCZ 
r n 

1 + At k 2 C n H 
Ch -f At hi 
1 + At k 2 C n H 
C n F + At k 2 C^C n H 
1 + A tk x C^ 

CtfF + At C? 4- knC^C^) 


The energy equation is solved explicitly. To start the 
chain reaction 0.00002 mole fraction of F was released 
into the flowfield initially. 

Using the above equations to calculate the reac- 
tion in a laminar premixed flame, we found that the 
time required for this reaction is the order of one mi- 
crosecond (Figs. 1 and 2). The turbulence time scale 
in the calculation, r = k/c , is the order of one sec- 
ond. Considering this large difference of time scales, 
it is impractical to use finite rate chemistry in this 
calculation. 


Fast reaction 

Since the chemical reaction is very fast compared 
to the flowfield development and the reverse reactions 
are negligible, we chose to use the following complete 
irreversible reaction mechanism. We assume that af- 
ter molecular mixing, complete reaction is achieved 
within each particle during the time interval of one 
marching step in the calculation. The reaction is cal- 
culated using the following equations. 

For > C Fi : 

eg 1 = o 

/■^n+1 /•'»n . o/^n 

° HF “ °//F + 

For Q 3 < Cp a : 


/•▼n+1 /^n 

Op a - 

/ orr *+l 

° HF — ° HF + Z{ ^H 7 
The energy equation can be written as 

i/, = -*}„ A Cr (OBZ) . 

where Mhf if the molecular weight of H F. 

Since the concentrations of reactants are low, the 
heat release is low, and the temperature rise is within 


j-n + 1 _ ‘j-n _j_ 

C p 

The error involved in this approximation is less than 

5%. 


Results and Discussions 
Code validation. 

Before applied to the hydrogen-fluorine diffusion 
flames, the computer algorithm is first validated us- 
ing experimental data for a heated turbulent jet. The 
non-reacting flow data serves as a check for the con- 
vection and molecular mixing process in the pdf solver 
as well as the k — c model in the N.-S. solver. 

Extensive experimental results for turbulent 
plane jet have been reported by many authors. Mea- 
surements for mean velocity field in a turbulent jet 
were first reported in the 1930’s [9]; turbulent shear 
stress measurements had been reported more re- 
cently[10, 1 1] . To determine the effect of turbulence 
on mixing, the temperature field of a heated turbu- 
lent jet had been studied by several authors. The tur- 
bulent jet has a slightly higher temperature than the 
ambient. Measurements of both the the mean tem- 
perature and the rms of the temperature fluctuations 
were given[12-17]. 

In the present study of this non-reacting flow 
case, the temperature field is treated as a conserved 
scalar and is simulated by the pdf of the tempera- 
ture; the velocity field and turbulent shear stress are 
obtained by solving the N.-S. equation and a k — c 
turbulence model. 

In the finite difference solution of the flowfield, 
41 grid points are used across half of the jet width. A 
symmetry boundary condition is used at the jet cen- 
terline. Fig. 3 shows the comparison of the present 
solution of the mean velocity field with the experimen- 
tal data, and Fig. 4 presents the numerical solution 
of the turbulent shear stress as compared to the ex- 
perimental data. Good agreements between numer- 
ical results and experimental data are observed for 
both the mean velocity field and the turbulent shear 
stress. The turbulent shear stress is calculated from 
— < u'v' >= u T du/dy , where i/ T = Cpk 2 /c. A good 
prediction of the turbulent shear stress ensures that 
the turbulent time scale, r = k/c, supplied to the 
Monte Carlo simulation is correct. 

For the Monte Carlo simulation of the temper- 
ature field, two sample sizes of 1000 and 1500 sam- 
ple particles per cell are used. The predicted mean 
temperature and the root mean square (rms) of the 
temperature variation are presented in Figs. 5 and 
6. The results show that both calculations produce 
fairly good comparisons with the experimental data, 
which means a sample size of 1000 particles per cell 
is large enough.’ 
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The above validation lends credibility to the N.S. 
solver as well as the Monte Carlo soiver developed in 
the present study. 

— Fn diffusion flames 

The flow conditions for the Hn — Fi diffusion 
flames are set according to an experiment performed 
by Hermanson and Dimotakis (1989). The flame con- 
sists two streams. The upper stream contains Nn and 
F 2 , the flow velocity is U\ = 22 m/s; the lower stream 
contains of iV 2 and of // 2 , with velocity f/ 2 = 8.8 m/s. 
In the present study, 6 cases involving various per- 
centage of H 2 and F 2 in the upper and lower stream 
are considered; the conditions are listed in Table 1, 
where A T a ^ is the difference between the adiabatic 
flame temperature and the free stream temperature. 


case 

No. 

equivalence 

ratio 

lower stream, 
mole fraction 
of H 7 

upper stream, 
mole fraction 
of F 7 

A 

1 

1 

0.02 

0.02 

186 

2 

1 

0.04 

0.04 

368 

3 

1 

0.06 

0.06 

554 

4 

1/4 

0.01 

0.04 

151 

5 

1/4 

0.02 

0.08 

302 

6 

1/4 

0.04 

0.16 

600 


Table 1. Initial conditions of the flames calcu- 
lated in the present study. 


Fig. 7 shows the calculated temperature rises due 
to combustion for cases 2 and 5 of Table 1 and the cor- 
responding experimental data. In the figure. 6t is the 
shear layer thickness determined by 1% of the tem- 
perature rise. AT is the actual temperature rise due 
to combustion. ( the two streams have the same tem- 
peratures initially,) and A T a y is the adiabatic flame 
temperature assuming complete reaction. The solu- 
tion for case 2 agrees fairly well with the experimental 
data for the same case, while the mean temperature 
is slightly over predicted for case 5. The computation 
shows the flame center shifts as a result of change in 
equivalence ratio, which is consistent with the exper- 
imental results. 

The rms of the temperature variance for cases 2 
and 5 are given in Fig. 8. The results show that in 
a diffusion flame, the temperature variance has two 
peaks, and the highest values do not coincide with 
the maximum values of the temperature distribution. 
The shifting of the flame due to the change of equiv- 
alence ratio can also be observed from this figure. 

The mass fractions of /f 2 , F 2 and HF for cases 
2 and 5 are presented in Figs. 9 and 10, respectively. 
One can see that although a complete reaction chem- 
istry model was used, with the pdf method, results 
simular to that of a finite rate computation are pro- 
duced, which is one of the many advantages of the 
pdf method. 

One experimentally established fact is that for 
the same equivalence ratio, with various mole con- 
centrations of fuel and oxidizer in the flowfield, the 


normalized temperature stays the same [4,5]. The 
present calculation confirmed this. Figs. 11 and 12 
are the mean temperatures and rms’s of the temper- 
ature variance for an equivalence ratio of one; three 
mole concentrations of fuel were considered. One can 
see that the curves coincide with each other. The 
same agreements are found for equivalence ratio 1/4 
(Figs. 13 and 14). 

The pdf distributions for flame temperature at 
the center and at the outer edge of the flame (for case 
2 of Table 1) are plotted in Fig. 15, where the x- 
axis denotes (T— < T >)/cr, with < T > = r being 
the mean temperature and a the rms. The pdf dis- 
tributions show that at the center of the flame, the 
temperature with the highest probability is not far 
from the adiabatic flame temperature, while at the 
outer edge of the flame, most of the time one would 
find a temperature close to that of the free stream 
temperature. Nonetheless, as a result of external and 
internal intermittency, low temperature fluid does ex- 
ist at the center and high temperature fluid the outer 
edge. Pdf distributions for various species concentra- 
tion can also be obtained from the solution, but will 
not be presented here. 

Concluding Remarks 

A Monte Carlo solution algorithm for the pdf 
evolution equation has been developed and success- 
fully combined with a finite difference flow soiver in 
the study of turbulent combustion. The algorithm 
was validated using turbulent mixing data from non- 
reacting flows. Turbulent diffusion flames of Hn- 
F 2 were computed using the pdf method, and good 
agreements between numerical solution and experi- 
mental data were observed. The computation iden- 
tified the change of equivalence ratio as the cause of 
flame shift, and demonstrated that similarity solu- 
tion exists for flames with the same equivalence ra- 
tio. The present work showed that a grid depen- 
dent Monte Carlo scheme can be readily applied to 
a general curvilinear coordinate system: therefore, it 
is suitable for realistic reacting flow computations. 
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time (10~ 6 sec) 

Figure 1. Temperature rise in a premixed lami- 
nar Hi-Fi flame. 



time (10 6 sec) 


Figure 2. Changes of mass fractions in a pre- 
mixed laminar Hj-F? flame. tf 2 , — - 

HF. 



Figure 3. Mean velocity distribution in a 2D tur- 
bulent jet. o Ref (9], — Present calculation. 



y/y.su 0 

Figure 4. Turbulent shear stress in a 2D turbu- 
lent jet. o Ref [10], A Ref [11], — Present calculation. 



Figure 5. Mean temperature in heated plane jet. 
A Ref [14], n Ref [17], o Ref [16] o Ref [15]. — pdf 
solution. 1500 particles per ceil; — - pdf solution. 
1000 particles per cell. 

0.40 | 

0.35 - 



0.0 0.5 1.0 1.5 2.0 


y/ysT 0 

Figure 6. RMS of temperature variance in heated 
plane jet. o Ref. (12], o Ref [17], A Ref (16), pdf 
solution. 1500 particles per cell; - - - pdf solution. 
1000 particles per cell. 
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Figure 7. Mean temperature rises in H 2 -F 2 
diffusion flames. Case 2 (refer to Table 1): o Ref [5], 
— present; Case 5: A Ref [5], present. 



Figure 8. RMS of Temperature variance in H 2 - 
A diffusion flames. — equivalence ratio 1 (case 2, 
Table 1). equivalence ratio 1/4 (case 5). 



Figure 9. Mass fraction distributions in an H 2 - 

F 2 flame, case 2, equivalence ratio 1. — H 2y A 

-- HF. 



Figure 10. Mass fraction distributions in an H 2 - 

F 2 flame, case 5, equivalence ratio 1/4. — tf 2 , 

F 2t --HF. 



Figure 11. Mean temperature rises in H 2 -F~. 
flames, equivalence ratio 1. Referring to Table 1: — 
case 1, o case 2, A case 3. 



Figure 12. RMS of temperature variance in H 2 - 
F 2 flames, equivalence ratio 1. Referring to Table 1: 
— case 1, o case 2, A case 3. 


85 


\0 



Figure 13. Mean temperature rises in H 2 -F 2 
flames, equivalence ratio 1/4. Referring to Table 1: 
— case 4, o case 5, A case 6. 



Figure 14. RMS of temperature variance in H?- 
F? flames, equivalence ratio 1/4. Referring to Table 
1: — case 4, o case 5, A case 6. 



( T -<T>)/a 


Figure 15. Probability density function distribu- 
tions in an H 2 -F 2 flame (case 2). — pdf at the center 
of the flame, pdf at the out edge of the flame. 
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An Improved k-e Model for Near Wall Turbulence 

T.H. Shih* and A.T. Hsu** 

Center for Modeling of Turbulence and Transition 
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Abstract 

This paper presents an improved k-e model for low Reynolds number turbulence 
near a wall. The work is twofold: In the first part, the near-wall asymptotic behavior of 
the eddy viscosity and the pressure transport term in the turbulent kinetic energy equation 
are analyzed. Based on these analyses, a modified eddy viscosity model with the correct 
near-wall behavior is suggested, and a model for the pressure transport term in the k- 
equation is proposed. In addition, a modeled dissipation rate equation is reformulated, 
and a boundary condition for the dissipation rate is suggested. In the second part of the 
work, one of the deficiencies of the existing k — e models, namely, the wall distance (e.g., 
y + ) dependency of the equations and the damping functions, is examined. An improved 
model that does not depend on any wall distance is introduced. Fully developed turbulent 
channel flows and turbulent boundary layers over a flat plate are studied as validations for 
the proposed new models. Numerical results obtained from the present and other previous 
k-e models are compared with data from direct numerical simulation. The results show 
that the present k-e model, with added robustness, performs as well as or better them other 
existing models in predicting the behavior of near-wall turbulence. 

1. Introduction 

The k-e model is one of the most widely used turbulence models in engineering 
applications. Patel et al.W recently reviewed existing two-equation models that can be 
integrated directly to the wall. One of their conclusions was that the damping functions 
used in turbulence models, especially the one for the eddy viscosity, need to be further 
modified in order to improve model performance. In fact, as we shall see later, many 
existing k-e models do not provide the correct near- wall behavior of the eddy viscosity. 

Shiht 2 ! recently proposed a new near- wall k — e model based on asymptotic analysis. 
The present paper is a direct extension of that work. 

In the present paper, we will first analyze, in section 2, the near-wall asymptotic 
behavior of the eddy viscosity and the pressure transport term in the fc-equation, and 
in sections 3 and 4, propose models according to their near-wall behaviors. The model 
equation for the dissipation rate is reformulated following an argument similar to that of 
Lumley,l 5 l and a boundary condition for e is suggested. 

An asymptotic analysis shows that, in the near wall region, while the pressure 
transport term in the turbulent kinetic energy equation is smedl compared to the dissipation 
and molecular diffusion terms, it is much larger than the turbulent transport term, and 
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** Supervisor, Comput. Phys. Section, Sverdrup Technology, Inc., Member AIAA. 

This paper is declared a work of the U.S. Government and 
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it plays an important role in the balance between the dissipation and molecular diffusion 
terms. This near-wall behavior is also observed in direct numerical simulation of fully 
developed channel flows (Mansour et alJ 32 ^, Kim et alJ 4 ^). However, in existing k-c models, 
this pressure transport term is either ignored or lumped into the turbulent transport model. 
The present work introduces a model for the pressure transport term explicitly. 

Most of the existing k-e models for near- wall turbulence use y + (defined as u r y/i/, 
where u T the friction velocity) as a parameter in constructing damping functions, with the 
Jones-Launder model being the only exception. While the use of y + is perfectly fine for 
simple attached boundary layer flows, it is inconvenient in more general applications such 
as separated flows and flows with corners, where y 4 * is not well defined. The Jones-Launder 
model has the advantage of avoiding y + ; however, the model is known to perform poorly 
in predicting near wall turbulent quantities, especially the turbulence kinetic energy. In 
the present work, a new damping function is derived based on asymptotic analysis (Section 
5). The new function is constructed upon a non-dimensional quantity that is independent 
of the coordinate system. 

The new models proposed in this paper were validated using direct numerical simu- 
lation data for fully developed turbulent channel flows and turbulent boundary layers over 
a flat plate. These numerical results are reported in Section 6. Comparisons are also made 
with other popular k — e models implemented in the same computer code. The numerical 
results show that the present model, in general, performs better than the existing models 
while providing added robustness. 

2. Asymptotic Analysis 

To analyze the near-wall asymptotic behavior of the eddy viscosity and other tur- 
bulent quantities, we expand the fluctuating velocities and pressure in Taylor series about 
the wall distance as follows: 


u \ = biy + C\y 2 + diy 3 + . . . 

= c 2 y 2 + d 2 y 3 + . . . 

u 3 — hy + c 3 y 2 + d 2 y 3 -t- . . . 

p = a p + b p y + c p y 2 + d p y 3 + . . . 


( 1 ) 


where the coefficients a p ,bi,c 2 , ... axe functions of x,z and t. Using the continuity and 
momentum equations, Mansour et al.^ derived the following relations between the coeffi- 
cients, 

2c 2 = — (&i,i + 63,3) 

dpi = 2isci (2) 

a Pt3 = 2uc 3 

The eddy viscosity is usually defined as 


(u,uj) — v T (U ltj + Uj t i) 


( 3 ) 
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whore ( ) stands for ensemble average and k = {v.a t )/2 is tlie turbulent kinetic energy. 
For plane shear flows, we can write from Eq. (3) 


— (uv) 
" r = dU/dy 


(4) 


and using Eq. (1), we obtain the near-wall asymptotic behavior of the eddy viscosity: 

vr = + dMi + C * C2 > + . 2 _ (^^ ), ( £ i | t/ 4 + 0(y (5) 

(pi) (M 

That is, near the wall vj is 0(y 3 ). A correct eddy viscosity model should have this near- 
wall behavior. We shall see later that many existing models do not have this near-wall 
behavior. For later use, let us examine also the near-wall asymptotic behavior of the 
turbulent kinetic energy k and its dissipation rate e = v ( u i,j u i,j)- Using Eq. (1), we 
obtain the following relations for the k and e: 

* “ ^ ^ 2 ^ + (^3^3 ))y 3 + 0 ( y 4 ) ( 6 ) 

= (&i) + (& 3 ) + 4((6 !Ci) + (& 3 c 3 ))y + 0(y 2 ) (7) 

V f 

In addition, the pressure transport term in the A:-equation, II = — -(uiP f i), becomes (using 
Eq.s (1) and (2)) 

n = -2u((b l c 1 ) + (b,c 3 ))y + 0(y 2 ) (S) 

The turbulent transport term in the fc-equation, — can be estimated as 0(y 3 ). 

Therefore, the pressure transport term is much larger than the turbulent transport term 
near the wall. 


3. Eddy viscosity model 

In this section, we will propose a model for the eddy viscosity using its near-wall 
behavior described in the previous section. In general, the eddy viscosity model can be 
written as 

v T = c v! (! (9) 

where v! and (! are the turbulent characteristic velocity and length scale, respectively. 
Depending on how the velocit}’ and length scales are specified, the eddy viscosity model 
can be a mixing length model, a one-equation ( k ) model or a two-equation (e.g. k-e) 
model. For example, in plane shear flows, Prandtl’s mixing length model specifies the 
characteristic velocity with t! dU/dy. For near wall turbulence, the Van Driest mixing 
length model further damps the length scale to y[l — exp(— y + /A)]. For more advanced 
mixing length models, see Baldwin and Lomax^, and Kingt 7 !. One-equation ( k ) models 
use k as the characteristic velocity, which is determined by the turbulent kinetic energy 
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equation. In two-equation models, e.g. k-e models, the length scale is usually specified 
by A: 3 / 2 / e, where e is determined by a dissipation rate equation. In this paper we will 
concentrate on two-equation models, wherein the eddy viscosity is usually written as 


k 2 

V T = 


( 10 ) 


where = 0.09, and f ^ is a damping function. The form of the damping function 

is critical in such formulations, since the prediction of the mean velocity field depends 
primarily on the eddy viscosity model. Thus it is important for an eddy viscosity model 
to have the proper near- wall behavior. We have examined the near- wall behavior of eddy 
viscosity models based on various k-e model equations. The results are listed in Table 1, 
which shows that some of the k-e models do not have the correct near- wall behavior of the 
eddy viscosity, namely, v t = 0(y 3 ). 

The quantity A: 3/2 /e is usually considered a characteristic length scale, , (or the 
size) of the energy containing eddies. One expects that near the wall the size of these 
eddies to be of the order of the distance from the wall, 0(y). However, Eq.s (6) and (7) 
show that k 3 / 2 /e is 0(y 3 ). Hence, k 3 ! 2 je is not an appropriate quantity to represent the 
length scale of the large eddies near the wall. We therefore introduce a new variable e: 


e 


e 


v 


dk/dx{ dk/dxi 

2k 


(ii) 


which has the following property: e approaches e away from the wall and is 0(y 2 ) near the 
wall. Therefore, fc 3 / 2 /e is a proper quantity to characterize the length scale of the large 
eddies. With this length scale, the eddy viscosity should be written as 


k 2 

U T = C^f^ — 


( 12 ) 


Now in order for v? to have the correct near- wall behavior, the damping function f ^ 
must be 0(y) near the wall and approaches 1 away from the wall. The damping functions 
used in various k-e models are listed in Table 2. If we consider the presence of the wall 
as the main effect on the eddy viscosity, then we may assume f ^ is mainly a function of 
y + . The form of f ^ can be determined quite accurately if we know i/^, k and e from, for 
example, the direct numerical simulation. One may optimize the following simple form by 
numerical experiments: 


f ti = 1 exp(-<x!y + - a 2 y +2 - a 3 y +3 - a 4 y +4 ) (13) 

The optimal values for channel flows axe a\ = 6 x 10 _3 ,a2 = 4 x 10 “ 4 ,a 3 = —2.5 x 
10 _6 ,a 4 = 4 x 10~ 9 . It can be shown that this form of damping function does provide 
the required near-wall behavior. As will be shown later, the above constants are valid for 
general boundary layer flows. 
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4. Modeled k-c equation 

To complete the eddy viscosity model, we need the modeled equations for the tur- 
bulent kinetic energy and its dissipation rate. In this section we will analyze the near-wall 
behavior of the Adequation and propose a model for the pressure transport term with a 
proper near-wall behavior. The equation for the dissipation rate is also reformulated with 
a formal invariant analysis. 


4.1 A:-equation 

We start with the equation for the turbulent kinetic energy, 

k,t + Ujkj = D y + T* + n + .P — t (14) 

where D v , T and II represent the transport of the turbulent kinetic energy due to the 
viscosity, turbulent velocity and pressure, respectively. P and e axe the rate of production 
and dissipation of the turbulent kinetic energy. The terms on the right hand side of Eq. 
(14) are defined a s follows: 

D v — vkjj 
T — —{kuj)j 


n = --(puj)j 

p 


(15) 


P = -(u,Uj)UiJ 
e = u{u itj uij) 

Using Eq.s (1) and (2), we obtain the budget of the Adequation near the wall, 

i=°(^ 

D„ = i/({6 2 ) + (b\)) + 6^((6 1 c 1 ) + (b 3 c 3 ))y + 0(y 2 ) 

T = °(y 3 ) (16) 

n = — 2i/((6iCi) + {b 3 c 3 ))y + 0(y 2 ) 

P = 0(y 3 ) 

e = u((b]) + (bl)) + M(6 lCl ) + (b 3 c 3 ))y + 0(y 2 ) 

This budget shows that the term II is much larger than the term T, and that II cannot 
be neglected if we want the Adequation be balanced in the near-wail region. However, the 
existing models either do not consider this term or simply combine it with the term T and 
model them as 


-( ku i)j = 


^ k 


(17) 


In this paper, we propose a model for IT which has a similar form to that of the standard 
turbulent transport model, but with a coefficient to ensure its correct near -wall behavior, 
Eq. (8). The proposed model form of n is 

Co 


n = 


l/T 


/^[l -exp( -?/+)] o k 


Ki 


(18) 
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where Co = 0.05 is a model constant. In some existing k-e models, it is assumed that e = 0 
at the wall. In that case, in order to balance the term D u , a nonzero artificial term D 
must be added to the ^-equation. The form of D for various k-e models is listed in Table 
3. Finally, the modeled fc-equation becomes 


dk rT dk d r/ irr.dk, „ ,dU { dU^dUi 

at + u ’W i ~ ar} (v + 71^} + n + + a ^- 'alj " ' + D 


(19) 


In the present model, D = 0, since e is nonzero at the wall. (The boundary condition for 
e will be discussed later.) 

4.2 e-equation. 

The exact dissipation rate equation is 


€ , + Ujej = D € u + T € + U e + PD (20) 

where D e u , T € and IT represent the diffusion rate of the dissipation rate due to the viscosity, 
turbulent velocity and pressure, respectively, and PD stands for the entire mechanism of 
the production and destruction of the dissipation rate e. The terms on the right hand 
side of the above equation are as follows: 


D l = i,e jj 

T ( = -u{u itk u itk u } ) yj 

n< = ~( P,kUj, k ),j (21) 

PD — 2v{{u i >k Uj tk ) -h )C7j i , 2u(ujUi tk )Ui tk j 

2i/(u, i jtUj > jtu, iJ ) 2v {ui' k jUi' k j) 

The term II f is usually neglected and the term T ( is modeled as 


T‘ = 



( 22 ) 


To model PD, we define $ by 

„ ^ 66 
PD = — — 4> 

At the level of the k-e model, we assume ^ is a function of v , vy, k, e x e, Uij and Uijk- 
Since ^ is an invariant, it must be a function of the invariants that can be constructed 
from these quantities. Therefore we can write 


* = <HRu " TUl : iUi - i 

6 ' J ee 


where R t is the turbulent Reynolds number k 2 / ve. We expand ^ in a Taylor series about 
i 'rUijUij/e and v vr Ui j * Uij kk/ee, and take only the linear terms. We obtain 


$ = V> 0 + rpi + \1> 2 vvtU i j k Ui j k -z 

t ee 


(23) 
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where the coefficients xpo, t/’i and xp 2 arc in general functions of R t . Finally, the modeled 
dissipation rate equation becomes 


de de 

dt + ; dxj 


+ + - c 2 f 2 j + E (24) 


where, C\ and C 2 are the model constants, and f\ and /2 axe functions of R t . The term 
E in the present model comes from the last term in Eq. (23): 


E = vi'rUijkUijk (25) 

where we have taken xp 2 = — 1- The form of E and C 1 , C 2 , f\ and / 2 for various k-e models 
are listed in Tables 3 and 4. 


4.3 Boundary Condition for e. 

Many of the earlier k — e models use e = 0 as the boundary condition for the 
dissipation rate. It is now generally agreed that this is not the physically correct boundary 
condition. However, controversy still exists in what boundary condition should be used for 
the ‘dissipation rate. In some calculations, de/dy — 0 is used, which clearly has no physical 
background. Most models use the second derivative of the turbulent kinetic energy at the 
boundary as the boundar} r condition for e, as listed in Table 1. This condition comes 
directly from the fc-equation and is physically correct; however, it makes the problem very 
stiff and thus put very stringent restrictions on the choice of initial profiles for k and e. 
If the initial profile for k is not given correctly, the second derivative of k can become 
negative and cause the solution procedure to diverge. 

We propose the following boundary condition based on the asymptotic analysis. At 
y = 0, it is obvious from eqs. (6) and (7) that 


e = 2u 



(26) 


This expression is exact at the wall, and it does not add stiffness to the solution procedure. 

5. Deficiencies and Improvements of Existing Models 
5.1 Damping functions 

One of the deficiencies of the existing near wall models is that most of the wall 
damping coefficients are functions of a wall coordinate, such as the y + . This types of 
damping function works well only in the cases of attached boundary layer flows with simple 
geometries where y + is well defined. For practical engineering problems with corner flows 
or separated flows, some ad hoc treatment to the damping function must be made. The 
same is also true for the length scale in an algebraic model. The only exception to the above 
is the Jones-Launder model in which the damping function is a function of R t = k 2 /ue. 

Although the Jones-Launder model has the advantage of independent of y + , it 
is known that this model does not predict correctly the near-all turbulence, especially, 
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the near-wall turbulent kinetic energy is underpredicted. One natural option could be to 
modify the J-L model such that it would predict the correct near- wall turbulent quantities. 
However, one basic difficulty met in such attempts is that the -equation is very stiff. The 
fact that we are putting k 2 back into the equation by using the parameter R t aggravates 
the situation. 


To avoid the above difficulty and yet still achieve the same end, we introduce the 
following parameter: 


R u 



(27) 


where U is the total velocity. (Note: U is the total velocity in a coordinate frame fixed to 
the solid boundary.) From the results given in Section 3, since e approaches a finite value 
at the wall, it is obvious from eq. (1) that R u = 0(y 4 ) near the wall. Similar to eq. (13) 
of Section 3, we write the damping function for the eddy viscosity as: 


fr = 1 - exp{- ai R l J 4 - a 2 R l J 2 - a 3 R u ) (28) 

with ai = 5 X lO” 3 , a 2 = 7 x 10“ 5 , and a 3 = 8 x 10 -7 . One can easily verify that with this 
damping function, the eddy viscosity has the correct near wall behavior, i.e., v? = 0(y 3 ). 

One point worth mentioning is the wide applicability of the above damping func- 
tion. Though developed with Shih model (Section 3, 4) in mind, it can be used with any 
existing k — e model that uses a non- zero boundary condition for the dissipation. This new 
parameter R u , unlike i? t , by no means affect the stability of the solution procedure. 

5.2 Pressure transport term. 

In order to remove the coordinate dependency of the fc-equation, we replace the 
pressure transport term given in eq. (18) by the following expression: 


(29) 
J 

where the model constant is readjusted to Co = 0.01. 

5.3 The formulation of e. 

In order to obtain the correct wall behavior for the eddy viscosity, we have in- 
troduced e in eq. (11), Section 3. Theoretical analysis shows that e is always positive. 
However, in numerical calculations, the value of e may become negative or even oscillatory 
due to round of errors. ( Depending on the accuracy of the numerical procedure, this may 
or may not be the case.) Here, an alternative definition of l is suggested: 


n = 


Cp VT y 
SI o k 


€ = (l-eip(-i? t 1/2 ))e 


(30) 


This expression has the same near wall behavior as eq. (11) but is less likely to cause 
numerical instability. 
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With the modifications suggested above, the model constants need slight adjust- 
ment. The constants used in the present modified model axe a s shown in Table 4. 

6. Numerical Testing 

Flows with self-similar solutions are particularly useful for accurate model testing, 
because their solutions are independent of initial conditions, and one does not need to 
choose carefully the initial conditions for the k and e. In this paper, we use a fully 
developed channel flow and a flat plate boundary layer for model testing. These flows are 
the simplest wall boimded turbulent shear flows with self-similar solutions. However, the 
complex features of the turbulence, for example, the effect of the wall on shear turbulence, 
are present. In the case of the channel flow, the k-e equations form a one-dimensional 
problem, numerical calculations are easy and accurate. Recently, the measurements^ 8 ' 
confirmed the accuracy of the direct numerical simulation dataJ 4 ' These data are used for 
model validations. 

In legends of the figures presented at the end of the paper, which will be discussed 
in detail in the following paragraphs, the word ‘‘present’ 1 refers to results obtained using 
the model suggested in Section 4 together with the new damping function and 6 given in 
Section 5, while the label “Shih” refers to results obtained using strictly formulations given 
in Section 4. 

6.1 Fully developed turbulent channel flow 

Let h be the half width of the channel, u T the friction velocity and Re r the Reynolds 
number defined as u T h/v. Let J7, k , 6, v x and y be the non-dimensional quantities, normal- 
ized by u r ,u^, u\/h,v and /i, respectively. The modeled equations for the channel flow 
become 


^=Re r 

dy 1 + i' 


d t 1 r . i dk \ t dU \ 2 1 

+ + rT t 


ux^dk 

i j ~r 

Ok dy 


-e = 0 


fdU. 2 1 „ , 66 ,<PU. 2 


1 


d f 1 r'T v dc 6 f klw v i x _ cc / u. \ j 

dy^Re^ + + HI~ T ~ 2 * 2 ~k + Rt\ 


(31) 

(32) 
0 (33) 


where 


k 2 

vt — C n fn Rc r ~z~ 

- = 

6 e 2 A: Re T 

fn = equation(13) 

/2 = l- Ii e X P [-( — 

c = C ° 

/m( 1 -exp(-y+)] 


) 2 ] 


(34) 
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The boundary conditions are simple. At. ihe wall, 

U = k = 0 

(g ) 2 

2k Re T 

and at the center of the channel, 

dk de 
dy dy 


( 35 ) 


(36) 


The main results from different k-e models for a channel flow with Re r = 180 
are plotted in figures 1-4. All the calculations are compared with the direct numerical 
simulation data. Figure 1 shows mean velocity profile, Figure 2 shows the turbulent kinetic 
energy, Figure 3 shows the turbulent shear stress distribution, and Figure 4 shows the 
Dissipation rate. From these numerical results we reach the following conclusions: The 
model of Jones and Launder^ 10 ) (JL) underpredicts the mean velocity as well as the peak 
value of the turbulent kinetic energy. Chien’s model! 1 d performs better than the JL model, 
but it overpredicts the mean velocity near the center of the channel as well as the turbulent 
kinetic energy. In these two models, e = 0 at the wall is used as the boundary condition, 
so the dissipation rate'near the wall cannot be correctly predicted. Lam and Bremhorst! 12 ! 
use a nonzero boundary condition for t and have made some improvement for the mean 
velocity and turbulent kinetic energy compared with the results of the JL model. However, 
the shear is much overpredicted, and the dissipation rate near the wall is not correct. The 
model of Nagano and Hishida! 13 ! presents a very good prediction for the mean velocity and 
shear stress, while the peak value of k is underpredicted. Their main modification to the 
JL model is a change in the damping function / M and the form of E. A zero dissipation rate 
at the wall is used. The numerical results from the Shih model and the present modified 
Shih model show improvements in the prediction of all quantities, including the dissipation 
rate. 

6.2 Boundary layer flows 

The boundary layer equations and the corresponding k — and e— equations are solved 
using a conventional semi-implicit finite difference scheme. In this scheme, the coefficients 
for the convection terms are lagged one step in the x-direction, and the source terms in 
the k— and e — equation are linearized in such a way that stability is ensured. 

In the present study, a grid of 100 points in the y-direction is used. The grid is 
stretched linearly with Ayj+\/ Ayj = 1.05. The grids expands in the y-direction according 
the the boundary layer growth rate. 

The results of this calculation are presented in Figure 5 through Figure 10. 

Figure 5 and 6 show the comparison of the wall shear stress from various models to 
experimental data and some DNS data. As shown in Figure 5, at low Reynolds number 
(based on momentum thickness), J-L model overpredicts the shear stress while Chien 
model and NH model underpredict the shear stress. One common character of these three 
models is that they all used zero boundary condition for the dissipation, and thus unable 
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to predict the correct turbulence near wall behavior. At high Reynolds number (Figure 
6) the JL model still overpredicts the wall shear stress, while the others seem to do a fare 
job. 

Figure 7 and 8 are the results for Reg = 1410, and Figure 9 and 10 are the results 
for Ree = 7700. From these results one can reach similar conclusions as we did from the 
channel flow case: The JL model generally underpredicts the peak value of the turbulent 
kinetic energy while overpredicts it in the inertia layer; the model also underpredicts the 
mean velocity profiles. The LB model and the NH model also underpredict the peak value 
of the turbulent kinetic energy near the wall. The three models mentioned above either 
have zero boundary condition for e or do not have the correct order of magnitude for eddy 
viscosity. The Figures show that, in general, the present model performs better than the 
existing models. 

Conclusions 

From the model testing, we conclude that the present k -e model has made consid- 
erable improvement over previous k-e models according to the comparison with the direct 
numerical simulation data. We find that the improvement is mainly due to the modified 
eddy viscosit}' model and the model of the pressure transport term in tho /c-cquation. The 
proposed dissipation rate equation also shows a better near- wall behavior than the previ- 
ous ones. The correct boundary condition for e also seem to play an important role in the 
accurate prediction of the turbulence near wall behavior. 
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Table 1 Eddy viscosity and boundary condition for e in various k-e models 



Model 

V r p 

BC: e w 


JL 

0(y 3 ) 

0 


Reynolds 

o(y 5 ) 

u dUi 

V dy* 


LB 

0(y 4 ) 

u §Us. 

U dy* 


Chien 

0(y 3 ) 

0 


NH 

0(y 4 ) 

0 


Shih 

0(y 3 ) 



Present 

0(y 3 ) 



Table 2 Damping functions used in various 

k-e models 

Model 


U h 

fi 

JL 


ex P ( \ + R t /SO ) 1-0 

1 — .3exp(— R 2 ) 

Reynolds 

1 - exp(-.019Si? t ) 1.0 

(1 - ,3exp(-ii 2 /9)]/^ 

LB 


[1 - exp(-.0l65Jl*)] 2 l + (^) 3 

x(l + ¥) ' 

1 - exp( — Jijf ) 

Chien 


1 — exp(— .01 15y + ) 1.0 

1 — ,22exp(— i? 2 /36) 

NH 


[1 — exp(— y + /26.5)] 2 1.0 

1 — .3 exp (— R 2 ) 

Shih 


Eq. (13) 1.0 

1 - ,22exp(-R 2 /36) 

Present 


Eq. (28) 1.0 

1 - .22exp(— R 2 /36) 


Table 3 Model terms in various k-e models 


Model 

n 



D 

E 



JL 

0 



M S £Y 

to 
^ \y 

¥) : 

2 

Reynolds 

0 



0 

0 



LB 

0 



0 

0 



Chien 

0 



2i/K 

y 2 

-^exp(- 

•5y + ) 

NH 

0 



-M^Y 


u 

)(0) 2 

Shih 

Eq. 

(18) 


0 

-r(0 

Y 


Present 

Eq. 

(29) 


0 

-r(0 

Y 



Table 4 

Model constants in various 

k-e models 



Model 

c h 



C 2 

<? k 



JL 

.09 


1.45 

2.0 

1.0 


1.3 

Reynolds 

.084 


1.0 

1.83 

1.69 


1.3 

LB 

.09 


1.44 

1.92 

1.0 


1.3 

Chien 

.09 


1.35 

1.8 

1.0 


1.3 

NH 

.09 


1.45 

1.9 

1.0 


1.3 

Shih 

.09 


1.45 

2.0 

1.3 


1.3 

Present 

.09 


1.5 

2.0 

1.3 


1.3 

R t = I< 2 /ue, R k 

= y/Ky/v, 

y+ = 

u T yjv. 






99 


20 


+ 


15 


10 - 


/ 

r / 

— — Jones &c launder 

— — — Lam & Bremhorst 

Chien 

Nagano & Hishida 

Shih 

Present 

O O O DNS Data 

1 1 1 1 — i — i — i — i — i i i i i i i i i i 

10 1 10 Z 

y + 

Figure 1. Mean velocity profile for channel flow, Re r = 180 
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Figure 3. Turbulent shear stress for channel flow, Rt T = 180. 
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Figure 4. Dissipation rate for channel flow, Rt r = 180. 
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Figure 5. Wall shear stress coefficient, flat plate boundary layer flows, low Reynolds 
number data. 
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Figure 6. Wall shear stress coefficient, flat plate boundary layer flows, high Reynolds 
number data. 
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Figure 7. Turbulent kinetic energy for flat plate boundary layer, Reg = 1410. 




Figure 8. Mean velocity profile for flat plate boundary layer, Res = 1410. 
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Figure 9. Turbulent kinetic energy for flat plate boundary layer, Reg = 7700. 
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Figure 10. Mean velocity profile for flat plate boundary layer, He® = 7700. 
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Low Reynolds Number Two-Equation 
Modeling of Turbulent Flows 


V. Michelassi* T.-H. Shih 
Center for Modeling of Turbulence and Transition 
NASA Lewis Research Center 

March 28, 1991 


Abstract 

A new k — e turbulence model that accounts for viscous and wall effects is 
presented. The proposed formulation does not contain the local wall distance 
thereby making very simple the application to complex geometries. The for- 
mulation is based on an existing Ic — c model that proved to fit very well with 
the results of direct numerical simulation. The new form is compared with nine 
different two-equation models and with direct numerical simulation for a fully 
developed channel flow at Re = 3300. The simple flow configuration allows a 
comparison free from numerical inaccuracies. The computed results prove that 
few of the considered forms exhibit a satisfactory agreement with the channel 
flow data. The new model shows an improvement with respect to the existing 
formulations. 
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N.N. Man sour 
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T.H. Shih* 
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Lewis Research Center 
Cleveland, Ohio 44135 


ABSTRACT 

This paper presents an improved k-e model and a second order closure model for low- 
Reynolds number turbulence near a wall. For the k-t model, a modified form of the eddy 
viscosity having correct asymptotic near- wall behavior is suggested, and a model for the pres- 
sure diffusion term in the turbulent kinetic energy equation is proposed. For the second order 
closure model, we modify the existing models for the Reynolds-stress equations to have proper 
near-wall behavior. A dissipation rate equation for the turbulent kinetic energy is also refor- 
mulated. The proposed models satisfy realizability and will not produce unphysical behavior. 
Fully developed channel flows are used for model testing. The calculations are compared with 
direct numerical simulations. It is shown that the present models, both the k-e model and the 
second order closure model, perform well in predicting the behavior of the near wall turbulence. 
Significant improvements over previous models are obtained. 
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ADVANCES IN MODELING 
THE PRESSURE CORRELATION TERMS 
IN THE SECOND MOMENT EQUATIONS 

Tsan-Hsing Shih and Aamir Shabbir 
Institute for Computational Mechanics in Propulsion 
and Center for Modeling of Turbulence and Transition 
Lewis Research Center 
Cleveland, Ohio 44135 

and 

John L. Lumley 
Cornell University 
Ithaca, New York 14853 


ABSTRACT 

In developing turbulence models, different authors have proposed various model con- 
straints in am attempt to make the model equations more general (or universal). The most 
recent of these are the realizability principle (Lumley 1978, Schumann 1977), the linearity 
principle (Pope 1983), the rapid distortion theory (Reynolds 1987) and the material indif- 
ference principle (Speziale 1983). In this paper we will discuss several issues concerning 
these principles and will pay special attention to the realizability principle raised by Lum- 
ley (1978). Realizability (defined as the requirement of non-negative energy and Schwarz’ 
inequality between any fluctuating quantities) is the basic physical and mathematical prin- 
ciple that any modeled equation should obey. Hence, it is the most universal, important 
and also the minimal requirement for a model equation to prevent it from producing un- 
physical results. In this paper we will describe in detail the principle of realizability, derive 
the realizability conditions for various turbulence models, and propose the model forms 
for the pressure correlation terms in the second moment equations. Detailed comparisons 
of various turbulence models (Launder et al. 1975, Craft et al. 1989, Zeman and Lumley 
1976, Shih and Lumley 1985 and one proposed here) with experiments and direct numeri- 
cal simulations will be presented. As a special case of turbulence, we will also discuss the 
two-dimensional two-component turbulence modeling. 
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ABSTRACT 

Some new developments in two-equation models and second order closure models will 
be presented. Two-equation models (e.g., k-e model) have been widely used in CFD for 
engineering problems. Most of low- Reynolds number two-equation models contain some 
wall-distance damping functions to acount for the effect of wall on turbulence. However, 
this often causes the confusions and difficulties in computing flows with complex geometry 
and also needs an ad hoc treatment near the separation and reattachment points. In 
this paper, a set of modified two-equation models is proposed to remove abovementioned 
shortcomings. The calculations using various two-equation models are compared with 
direct numerical simulations of channel flows and flat boundary layers. 

Development of second order closure model will be also discussed with emphasis on the 
modeling of pressure related correlation terms and dissipation rates in the second moment 
equations. All the existing models poorly predict the normal stresses near the wall and fail 
to predict the 3 dimensional effect of mean flow on the turbulence (e.g., decrease in the 
shear stress caused by the cross flow in the boundary layer). The newly developed second 
order near- wall turbulence model to be described in this paper is capable of capturing the 
near-wall behavior of turbulence as well as the effect of three dimension mean flow on the 
turbulence. 
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1. k-e model 

The two-equation model, especially k-e model, is still the most widely used model 
for computing engineering flows. We first list some of the commonly used two-equation 
models,! [ 1 1i[ 2 1i[ 3 1»[ 4 ]»[ 5 1»[ 6 ] and their predictions on the fully developed channel flows and 
boundary layer flows compared with the corresponding direct numerical simulations. Then, 
we propose a modified k-e model which does not contain any wall distance. The proposed 
k-e model has been also tested using direct numerical simulation data. 

The eddy viscosity v t is assumed in two-equation models as follows: 


k 2 

v T = C,f^ 

e 


or 

vt — Cpf^kr 


where r = k 2 fe 

The general k-e (or k-r) model equations are of the following forms: 


+ Ujkj 




+ II + vt (Uij + Uj t i)Uij — e + D 


J 


e ,t + Uj e ,j 
r ,< + U J T J 



+ (1 - C^u T (Uij + U jt i)Uij + ( C £ 2/2 - 1) 
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This table lists the model terms, damping functions and model constants appeared in 
various two-equation models 
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The following figures show the predictions on fully develop channel flow using various 
two-equation models compared with the direct numerical simulation dataJ 7 ! Ploted quan- 
tities include the mean velocity U, turbulent shear stress (uv), turbulent kinetic energy k 
and dissipation rate EPS e. The open circle represents direct numerical simulation, and 
the solid line represents the model prediction. 
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The figures below show the model predictions on flat plate boundary layer flow. A 
direct numerical simulation of boundary layer flow^ is used for comparison. The skin 
friction coefficient C f is also included in the comparisons. 



1 13 


Overall, Shih’s k-e model gives better prediction in both fully developed channel flows 
and flat plate boundary layer flows according to the comparisons with corresponding direct 
numerical simulations. However, this model has the same problem as the others, that is it 
contains the wall distance parameter y + , which is defined as 



where u T is the friction velocity. The difficulty would occure in some situations. For 
example, near the seperation point u T ap roaches zero and hence u t (through //x(y + )) will 
approach zero everywhere when this u r is used. Another example is the flow with complex 
geometry that the wall distance is not well defined. In the both cases, the ad hoc treatment 
is needed in the model implimentation. We notice that Jonse-Launder’s model^ 1 ) does not 
contain y + . However, its present form does not perform very well in the simple testing 
flows. Here, we based on the Shih’s k-e model modify the parts which is related the y + 
with another parameter R = R is a ratio of turbulence length scale to viscous 

length scale. The modifications^ 9 ! made here are: 

u 

n 

€ 

where R t — k 2 /ve, Cq — .004, C 3 = .0004, C$ — 1.2 


= 1 — exp {C 3 [l-expCCefl 1 / 4 )]} 


——k 


= € 


in* "I,- 

1 — exp(— 
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turbulent viscosity 


The following figures show the predictions from the present model on fully developed 
channel flow compared with other models (including Jonse and Launder’s model^) and di- 
rect numerical simulation dataJ 7 ! The open symbols represent direct numerical simulation, 
and the lines represent the model prediction. 
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2. Second order modeling of near-wall turbulence 

Using the near-wall asymptotic behavior of turbulence^ 10 ! a 5 model constraints, we 
formed a set of modeled transport equations for the Reynolds-stress tensor and the dissi- 
pation rate of turbulent kinetic energy. The main emphasis was on developing a near-wall 
model for the pressure correlation and dissipation terms in the Reynolds-stress equation. 
A modeled dissipation rate equation is derived more rationally. Asymptotic analysis shows 
that near the wall, the viscous diffusion term in the Reynolds-stress equations becomes 
the leading term and is balanced by the pressure correlation and dissipation terms. We 
use this as a model constraint in the model development. The proposed models satisfy 
realizibility and ensure no unphysical behavior will occur. Here, we briefly describe and 
list the proposed models. 

Reynolds stress equation 

The exact equation for the Reynolds stress tensor is: 

§- t { «.'<*>> = Pn + Tn + + n y - 

where ( ) stands for an ensemble average, D/ Dt = d/dt + U k d/dx k . the terms P t j , T t; , 
d\j \ Iljj and £ij represent the production, turbulent diffusion, viscous diffusion, velocity 
pressure-gradient correlation and dissipation tensor, and are identified as follows: 

Pij = ~{uiU k )Uj tk ~ (ujU k )U it k 
Tij = -(uiUjU k ),k 
D if = V { U i U i),kk 
n u = -~ p ( u iP,i+ u iP,i) 

eij = 2is(u itk u jt k) 

The proposed near- wall model for II ij — €{j is: 

n,j - £ij = -f w j~[2(uiUj) + 4 ((uiUk)njTik + (u J u Jt )n,n J t) + 2(ittu,)njtn,n,nj] 

(?) 

where n, is a unit vector normal to the surface, and f w = exp(— (i?t/Ci) 2 ), Rt = 

C 1 = 1.35SR°J 4 , R er = u T S/v. u r is the friction velocity, S is the thickness of the boundary 
layer or the half width of the channel. 

Away from the wall, the velocity pressure-gradient correlation 11,-y is split into the 
rapid part 11^ and the slow part 11^: 

n« = "!>’ + njj> 
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The proposed model for II - ^ — e, ; - is: 

n\V-e ij = -^b ij + ^S ij )(l-f w ) 

where 

0 = 2 + + 8(U ln t X + 62A (~ n + 2.3///)]} exp( — ~ytz) 

y R t ' R/ 2 

F = 1 + 27/// 4- 9// 

II = - jM;i 
III = 

bij = ( UiUj)/(q 2 ) - <5 0 /3 

The rapid part of velocity pressure- gradient, II-^ is modeled as follows(Shih and 
Lumleyl 11,12 !): 


ng* = ( J + + u'i.i ) - |(1 - a s )(P,j - | PSi,) 

+ (1 + - \ ps a) + - Dii) + f M* 


3 3 

2 


[(( U » U *)^.M A ( u j u k)Ui t q)(v.kUq ) (tliU p ) (ujUg)(U p< g + t^ 9 l p)] 


where, 


/’ij = ~{uiV-k)Uj t k - ( UjU k )U it k 
Rij = {njHk)Uk t i 

P ~ ~P ■ 

2"* 

<■» = -- 5(1 + C 2 F‘^) 

C 2 = 0.8[1 - exp(— (i? t /40) 2 )] 

Finally the model for the third moments is modeled as: 

(uiUjUk) = -,07-|-^ [(u fc Up)(u i U J ) i p + (UjU p ){uiU k ) tP + (u.Up^tijU*)^] 
Dissipation rate equation 

The modeled dissipation rate equation derived in this work is: 


e,t + Uic ti = (t/e,,- - (eu,)) )t - - t/> 0 


€6 

(?y 


- - UwPijt 


1 18 


where 


xj> 0 = — + 0.98(1 - 0.33 ln(l - 55//)] exp(-2.83.R, 1/2 ) 
5 

ipi m 2.1 

2 = -.15(1- F) 

4(g a ) 

The turbulent flux term (euk) is modeled as: 

(eu k ) = -.07^-{u k u p )e iP 


These figures show some existing Reynolds stress models (for example, Launder and 
Shima^, Lai and So^ 14 ^) and present model compared with the direct numerical simulations. ^ 
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3. Second order modeling of a three-dimensional boundary layer 

A study* of three-dimensional effects on turbulent boundary layer were achieved by 
direct numerical simulation of a fully developed turbulent channel flow subjected to trans- 
verse pressure gradient. The results show that, in agreement with experimental data^ 16 !, the 
Reynolds stresses axe reduced with increasing three-dimensionality and that, near the wall, 
a lag develops between the stress and the strain rate. To model these three-dimensional 
effects on the turbulence, we have tried various two equation models and second order 
closure models. None of the current models can predict the reductions in the shear stress 
observed using direct numerical simulations. However, we found that the newly proposed 
second order closure model listed in the previous section do at least qualitatively capture 
these three-dimensional effects. 
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The following figures show the direct numerical simulation of the three dimensional 
boundary layer flow and the model predictions from Launder and Shima, Lai and So and 


the present models. 



Schematic of three-dimensional channel flow. 
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Turbulent intensities in fixed coordinates at A t 
0.15 intervals. 



Reynolds stress component in fixed coordinate at 
At = 0.15 intervals. 
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Distribution across the channel of the mean flow 

angle ( ), the Reynolds shear stress angle ( ) , 

and the mean velocity gradient angle ( ). 
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Abstract 

This paper presents a set of realizable second order models for boundary free turbulent 
flows. The constraints on second order models based on the realizability principle are re- 
examined. The rapid terms in the pressure correlations for both the Reynolds stress and 
the passive scalar flux equations are constructed to exactly satisfy the joint realizability. 
All other model terms (retum-to-isotropy, third moments and terms in the dissipation 
equations) already satisfy realizability (Lumley 1978, Shih and Lumley 1986). To correct 
the spreading rate of the axisymmetric jet, an extra term is added to the dissipation 
equation which accounts for the effect of mean vortex stretching on dissipation. The test 
flows used in this study are the mixing shear layer, plane jet, axisymmetric jet and plane 
wake. The numerical solutions show that the new unified model equations (with unchanged 
model constants) predict all these flows reasonably as the results compare well with the 
measurements. We expect that these model equations would be suitable for more complex 
and critical flows. 
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Evaluation of Turbulence Models 
for Predicting Buoyant Flows 

Experimental data for the buoyant axisymmetric plume are used to validate certain 
closure hypotheses employed in turbulence model equations for calculating buoyant 
flows . Closure formulations for the turbulent transport of momentum , thermal 
energy , kinetic energy , and squared temperature used in the k-t and algebraic stress 
models are investigated. Experimental data for the mean velocity, mean temperature , 
and kinetic energy are used in the closure formulation to obtain Reynolds stresses, 
heat fluxes, etc., which are then compared with their measured values. 


1 Introduction 

Various turbulence models have been formulated for pre- 
dicting buoyancy-driven flows. Some of the parameters in these 
models have been determined by keying the solution of the 
model equations to experimental data for certain basic flows 
such as decay of grid turbulence. Other parameters have been 
determined by calibrating closure formulations directly with 
experimental data. However, this approach may be somewhat 
inaccurate due to the lack of quality experimental data for 
certain correlations, especially dissipation. Finally, certain 
model parameters have been fine tuned or determined by re- 
quiring that the computed solution agree with experimental 
data for more complex flows, such as shear flows, in addition 
there have been instances where ‘model parameters have been 
adjusted or empirical corrective terms added so that agreement 
with experimental data is accomplished for a particular flow. 
When model parameters are adjusted to get agreement, say 
for the mean velocity and temperature fields for a particular 
flow, little regard is given for the internal integrity of the model. 
In other words, are the various processes such as diffusional 
transport, pressure-strain interactions, etc., predicted cor- 
rectly? Or are there compensating assumptions where one pro- 
cess is overpredictcd at the expense of another and yet the end 
predicted result for the mean flow agrees with experiment? 
The lack of complete sets of data for higher moments, dissi- 
pation, and pressure-velocity correlations for various flows 
has prevented detailed verification of closure models for the 
various processes that have to be modeled. 

The objective of this paper is to use the recently obtained 
and comprehensive experimental data of Shabbir and George 
(1987) and Shabbir (1987) on the axisymmetric buoyant plume 
to assess the various closure relations proposed for the kinetic- 
energy/dissipation and the algebraic stress models for buoy- 
ancy-dominated flows. The usual approach is to solve the 
modeled differential equations numerically, and then compare 
the computations with the experiment. However, this method 
docs not help pinpoint the drawbacks in the various terms of 
the models. In this paper, instead of the usual approach, cor- 
relations obtained from measured velocity and temperature 
arc used directly to verify the closure hypotheses for the tur- 
bulent transport of momentum, thermal energy, and turbulent 
kinetic energy. 

2 Experimental Data 

The data used were taken in an axisymmetric buoyant plume 
by Shabbir and George (1987) and Shabbir (1987). who meas- 
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ured velocity and temperature fields at several vertical levels 
above a heated source of air. Here we briefly summarize their 
experimental technique and results. 

The three-wire probe used consisted of a cross-wire and a 
temperature wire. Thus the instantaneous values of the two 
velocity components (vertical and radial) and temperature were 
measured. The axisymmeiry of the flow was established by 
using an array of 16 thermocouples and also by rotating the 
cross-wire by 90 deg. Profiles for the correlations between the 
velocity components and velocity components with tempera- 
ture through the fourth order were determined from the in- 
stantaneous measurements. 

Source conditions were continuously monitored in order to 
calculate the rate at which buoyancy was added at the source. 
The source Grashof number was 5.5. By integrating the mean 
energy equation, an integral constraint can be obtained for a 
buoyant plume. For a neutral environment this constraint im- 
plies that the rate at which buoyancy crosses each horizontal 
section is constant and must equal the rate at which buoyancy 
is added at the source, i.e., the ratio 

= jr |^2r g0(UAT+ut)rdr j (1) 

must be unity ( F c is the source buoyancy). This integral con- 
straint was satisfied within 7 percent. 

The correlation profiles at various heights were found to be 
similar in the coordinate tj = r/ 7 . ( z accounted for the virtual 
origin) when the velocity is scaled by U f Ff'z" 1 '* and the 
temperature is scaled by T s = Fy 2 z~ in /gp. The measurements 
agreed well with the earlier study by George et al. (1977), who 
measured only the temperature and the vertical component of 
velocity. The scatter in the measurements of higher moments 
is typical for such flows and is also present in previous ex- 
periments, such as those of George et al. (1977). The primary 
reason for the scatter is that slow time scales of the flow require 
much longer averaging time for the higher moments in order 
to obtain the same statistical convergence as for the mean 
quantities. Other errors in the measurements arise from the 
flow reversal on the hot wire — a phenomenon most likely to 
occur toward the outer edges of the flow where local turbulent 
intensities arc considerably higher. These arc discussed in Shab- 
bir and George (1987). 

The various correlations in similarity variables were fitted 
with curves using a least-squares fitting procedure. This rep- 
resentation allows easy evaluation of the terms in the governing 
equations and closure formulations when they arc cast in sim- 
ilarity variables. Using these profiles the balances for the mean 
momentum and energy differential equations were carried out 
to check whether the flow satisfied the equations of motion it 
is supposed to represent. Within the thin shear layer and the 
Boussincsq assumption the mean momentum and energy cqua- 
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Fig. 1 (a) Balances of mean energy and momentum equations (taken 
from Shabbir and George, 1987) 


tions can be respectively written as 

dU dU id 
U—+ V— =- -- (nw)-gPAT 
dz dr r dr 

rr dA T BAT 

u — — + v — — — (rvl) 

dz dr r dr 


( 2 ) 

(3) 


Since all the quantities appearing in these equations are meas- 
ured, their profiles were substituted to see whether the meas- 
urements balance the equations. Figure 1(d), taken from 
Shabbir and George (1987), shows that the experiment satisfies 
this nontrivial test within 10 percent. An error of such mag- 
nitude is typical of turbulent shear flows. 

The dissipation of mechanical energy was determined by 
balancing the turbulent energy equation 



Fig. 1(b) Mechanical and thermal dissipation profiles! Full lines are 
experimental; the chained line Is from model equation (10), which was 
solved for < with all other quantities taken from experiment. 





+ P+G-C 


(4) 


where P = — u-uj d(J/dxj is the mechanical production and 
G = is the production due to buoyancy. Each term 

except for the dissipation e and the pressure transport 'pu i is 
determined from the experimentally determined correlat ions . 
The pressure transport was evaluated from JnTj/ p = —q 2 ^/ 
5, a formula given by Lumley (1978). Although this closure 
relation has not been verified experimentally, it was felt that 
since the pressure transport is significant, some correction 
should be included rather than simply neglecting it, as is often 
done. The dissipation determined from the balance of the 
turbulent kinetic energy equation is shown in Fig. 1(a) as a 
solid line. 

By a similar procedure the dissipation of the mean-square 
temperature i 2 is determined from 


Ujj- = u/ —2ujl — 2t, 

1 dxj 3xj r 1 3xj 


(5) 


All terms are evaluated from experimental data and the re- 
sulting thermal dissipation is shown in Fig. 1(6). 

The time scales cf/t and r 2 / e , for the relaxation of the me- 
chanical and thermal dissipation, respectively, are shown in 
Fig. 2, along with their ratio 

R =(?/«,)/(?/< ) ( 6 ) 


Nomenclature 


F c = buoyancy flux, equation (I) 
g = acceleration due to gravity 
G = turbulence production from 
buoyancy 

k = turbulent kinetic energy 
p = fluctuating pressure 
P = turbulence production by mean 
flow 

Pr r = turbulent Prandtl number 


r = radial coordinate 
R = time scale ratio, equation (4) 

( = fluctuating temperature 
T = mean temperature 
u = fluctuating axial velocity com- 
ponent 

U = mean axial velocity component 
v = fluctuating radial velocity com- 
ponent 


z = vertical coordinate 
0 = coefficient of thermal expan- 
sion 

e = dissipation of mechanical en 
ergy 

e, = dissipation of mean-square 
temperature 

= turbulent eddy viscosity 
P = mean density 
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Transactions of the ASME 



Fig. 2 Mechanical and thermal time scale ratios and variation of R as 
calculated from experimental data 
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Fig. 3 Shear stress 
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Fig. 4 Axial heal flux 
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Fig. 5 Radial heat flux 


By invoking the thin shear layer assumption for a buoyant 
plume, the above relations reduce to 


_ dU 


ar 

ul = - (^/Pr 7 ) — 


These time scales appear extensively throughout the model 
formulations and will be further discussed in the next sections. 


3 Assessment of Closure Hypotheses of k-e Model 

The form of the k-t model, which is considered to be the 
standard one, is that used by Launder and Spalding (1974). 
In this model the Reynolds stress is given by 



and the heat flux by 


v T ar 


( 8 ) 


where v T - C^/t and = 0.09 (see Launder and Spalding, 
1974). 


vl = - (^ r /Pr r ) — 

Taking Pr^ = 1.0, the right-hand sides of the above equations 
were evaluated experimentally. These arc compared with meas- 
ured values of uu, ul, and T7 in Figs. 3-5. The points are 
experimental values and the chain lines arc from the model. 

The modeled values of uu and ul compare reasonably with 
the experimental profiles except in the outer portion of the 
curves. On the other hand the modeled profile of vertical heat 
flux ul is much smaller than the experimental one. It is well 
known that the simple gradient models given by equations (7) 
and (8) with an isotropic eddy viscosity arc inadequate for 
determining streamwise turbulent momentum and heat fluxes. 
Usually these quantities do not influence the prediction for 
shear layers since only the radial fluxes are important in these 
flows. However, in the case of the buoyant plume the flux 
ul is a dominant production term in the turbulent kincticenergy 
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Fig- 6 Transport ot kinetic energy 


equation and its correct calculation is very important for ac- 
curate prediction of k. 

The diffusional transport in the kinetic energy equation (4) 
for a thin shear layer is modeled as 


1—5 1 _ dk 


( 9 ) 


Using pv/p = -vcf/S from Lumley (1978) gives vq 2 /! = 
— (5/3)^ 7 dk/dr from which the result, with the right side eval- 
uated from experimental results; is showij in Fig. 6. It is seen 
that the predicted and experimental data peak at different 
radial locations; however, the predicted magnitude is more 
accurate, which indicates that the pressure diffusion needed 
to be taken into account. 

In the k-e model the dissipation is calculated from 

2 




( 10 ) 


where t'uj = - (vj/oJdt/dxj and o € = 1.3, C lX = 1.44, and 
C <3 = 1.92 as given by Launder and Spalding (1974). In order 
to get an indication of the validity of equation (10), it was 
numerically solved for the dissipation € with all other quantities 
needed to evaluate the coefficients determined from the ex- 
perimental correlations. The result is shown in Fig. 1(6) and 
it is seen that it compares reasonably well with the curve ob- 
tained from balancing the turbulent kinetic energy equation 
with experimental data. 

Launder et al. (1972) showed that the standard k-€ model 
yields a solution for the axisymmetric jet that overpredicts the 
spreading rate by about 30 percent. The standard k-t model 
also does not correctly predict the axisymmetric buoyant plume 
(Hossain and Rodi, 1982). Proposals have been made (Pope, 
1978; Hanjalic and Launder, 1980) for modifying the dissi- 
pation equation, based on arguments concerned with vortex 
or eddy structures characteristic of axisymmetric flows. The 
modified equation produces more dissipation, thus decreasing 
the turbulent eddy viscosity, which results in a smaller spread- 
ing rate of the flow. Here we use the empirical correction given 
by Rodi (1972) where C« 2 = 1.92 (1-0.035 H) with H - I (y E / 
U m )dU m /dx I 0 * 2 and where U m is the maximum velocity and 
y E is the distance from the centerline to the edge of the shear 
layer. This correction decreases the destruction term in the 
dissipation equation, hence producing an increased dissipation. 
However, when this correction is used in equation (10) there 
is very little change in the solution for < when experimental 
data are used for the other quantities in the equation. This is 
probably due to the approach taken here, which does not allow 
for the nonlinear interactions between the various terms in the 
closure. If the kinetic energy and dissipation equations arc 
solved simultaneously, then the axisymmetric correction will 
produce a significant change in the solution of the k-t model. 


4 Assessment of Closure Hypotheses for Algebraic 
Stress Model 

Chen and Rodi (1975), Tamanini (1978), Chen and Chen 
(1979), and Hossain and Rodi (1982) have made predictions 
for the buoyant jet using algebraic stress models. Many of the 
ideas used in these models for calculating buoyant flows orig- 
inated with Launder (1975, 1978). Algebraic stress models are 
obtained by simplifying the convective transport equations for 
Reynolds stresses and heat fluxes so they arc no longer dif- 
ferential equations. The dynamic equation for the Reynolds 
stress tensor is 

(C-D)— ^Pv+Gv- 




30C 2 — 22 


(ID 


where P tj - — u,u k d(J/dx k - uju k 3U/dx k is the mechanical 
production and G „ = — ujl — pgj is the buoyancy 
production. The left side represents convection minus diffu- 
sional transport, the dissipation is assumed isotropic, and the 
last three lines rep resent the closure formulation for 
p(du/dXj+ du/dxA/p given by Launder et al. (1975) and 
Launder (1975, 1978). Launder assumes: (1) an equilibrium 
situation where convection is balanced by diffusion (C - D 
= 0) and production is balanced by dissipation (P + G - < 
= 0); (2) the second and third terms (third line) in the rapid 
part of the pressure-velocity correlation are negligible; the 
coefficient C 2 is adjusted so that the first term approximates 
the entire rapid part; (3) the parameter C 3 is taken equal to 
C 2 . After applying all the assumptions 


u u - 

• j 


2 C, + C 2 - 1 
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kbjj 


1-Q 


2 k 


+ G,) (12) 


C, " C, 

where c, = 2.2 and C 2 = 0.6. It should be pointed out that 
in free shear flows the equilibrium condition (C — D = 0 and 
P + G - t = 0) only applies in the outer portion of the flow. 
Also, Zeman and Lumley (1976) found C 3 = 0.3, after applying 
all the constraints applicable to determining the contribution 
of buoyancy to the pressure-strain correlation. 

The dynamic equation for the heat flux is 


(C-D)- f = — Uftj ~~ -u? dU ‘ 


dXj 


_ < _ _ _,3G; 

C„-u7 + C 1< 5 7 — 


-Pgf 1 


1 _,dUj —j 

5 ^ BXi +Ci ^ 


(13) 


where the first line on the ri ght side is the production and the 
second line is the closure for pdt/dx/p. Neglecting convection 
and diffusion (C - D = 0) and the third term in the second 
line, equation (13) becomes 

(14) 

where C„ = 3.0, = 0.5, and C 3 , = 0.5. Zeman and Lumley 

(1976) show that = 0.8 and C 3/ = 0.2 from theoretical 
considerations. 

Neglecting the convection and diffusional transport in equa- 
tion (5) and eliminating with equation (6) gives 




(13) 
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which was given by Launder (1975, 1978) and used by Hossain 
and Rodi (1982). 

Chen and Rodi (1975) and Chen and Chen (1979) used the 
differential equation (5), with t t eliminated by using equation 
(6) to determine f 2 rather than using equation (15). In either 
ease R is taken to be a constant equal to 0.8 (Hossain and 
Rodi, 1982; Chen and Rodi, 1975; Chen and Chen. 1979; 
Launder, 1975, 1978). It is seen in Fig. 2 that the experimentally 
determined value of R is much lower with an average value 
across the profile of roughly 0.25. Launder (1978) cites ex- 
perimental evidence for R being in the range of 0.5 to 0.8. 
However, he found that the algebraic stress relations agreed 
best with an experiment for a stably stratified homogeneous 
shear flow with R = 0.8. Hence, that value has been adopted 
in the algebraic stress models. The experimental results of 
Shabbir and George (1987) indicate R is much lower for strongly 
buoyant flows. When the algebraic stress model is applied to 
this experiment with R = 0.8, the results arc very poor for 
r 2 . Therefore, in the following evaluation of the algebraic stress 
model, the experimentally determined profile for R (Fig. 2) is 
used. 

Equations (10), (12), and (13) represent a system of algebraic 
equations that can be solved for u,u Jt u 7, and t 2 . Employing 
the thin shear layer approximation, where only gradients in 
the radial direction are retained, Hossain and Rodi (1982) give 


^2 _ 2 C x + C 2 — 1 ^ + 1 — C 2 k 
“3 C, C, e 


/ dU \ 

f - 2 uv — + 2 figvt j 


(16) 



Fig. 7 Mean square temperature 
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(18) 
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from which v, = C,; tVt where 
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dT/dr\ 

dU/dr) 


(19) 

( 20 ) 

(20 

( 22 ) 


The system of equations (16)-(2I) was solved to determine the 
Reynolds stress and heat flux components with U, T, k , «, and 
R given by the experiment. The following values for the con- 
stants were used: 

C, = 2.2, C 2 = 0.6, C 3 =0.6 
C|,=3.0, C* = 0.5, C Jf = 0.5 

The value of C M , which appears in the eddy viscosity relation 
= CJfi/t and is given by equation (22), is roughly 0.125 
and is reasonably constant across the flow. This value is con- 
siderably larger than the value of C 0 = 0.09 in the standard 
k-t model. Rodi (1972), to correct for the discrepancies in the 
prediction for the axisymmetric jet, developed an empirical 
correction to C^. Thc parameter C„ is replaced by (1-0.465//) 
where H = 1 (yE/U m )dU m /dx\^ 1 ^U m is the maximum velocity, 
and y E is the distance from the centerline of the edge of the 
jet. Hossain and Rodi (1982), Chen and Rodi (1975), and Chen 
and Chen (1979) used this correction in their predictions for 
turbulent buoyant jets. When the correction is applied, we get 


approximately 0.09, which is the value of for the standard 
k-e model. 

The question we arc asking is, “given the turbulent energy, 
dissipation, velocity, and temperature, does the proposed al- 
gebraic stress expression correctly predict the Reynolds stress 
and heat flux components?** Figure 8 shows the radial Rey- 
nolds stress determined from equation (18). It is seen that the 
predicted value i^/k = 0.53 is a little smaller than the exper- 
imental curve in the center portion of the plume, but agrees 
quite well with experiment in the outer portion. The predicted 
shear stress Uv, the axial heat flux 177 , the radial heat flux F7, 
and the mean squared temperature r 2 arc shown in Figs. 3, 4, 
5, and. 7, respectively. Again the points are experimental data 
and the broken lines are from the model. It is seen that the 
shear stress uv and radial heat flux vl are predicted reasonably. 
However, the vertical heat flux 771 and temperature fluctuations 
F arc predicted poorly and have incorrect shapes; unlike the 
experimental values they go to zero near the origin. 

Equation (21) gives r 2 proportional to the radial temperature 
. gradient, which is zero at the centerline. Then since f 2 = 0 at 
r = 0, equation (19) gives w7 = 0 at r = 0. In order to obtain 
nonzero values, for ul and t 2 at. the centerline from the model 
equations (12), (14), and (15), terms containing the axial gra- 
dient, i.c., BU/dz and dT/dz t were retained. These terms were 
added to equations (19) and (21) and the system of equations 
was- solved again. Although the centerline values of ut and 
t 2 were found to be nonzero, the predictions still decreased to 
relatively small values near the centerline. 

Another possibility for this behavior is. the neglect of ad- 
vection and diffusion terms in the model. Gibson and Launder 
(1976) have proposed the following model for these terms: 


(C-D)^= ^ (P+C-t ) 


(23) 


(0-/7)^= ^ (i’,- < ,)+ ^ (P + C-t) (24) 

where P, is the production term in the l 1 equation. These were 
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incorporated in equations (16)— (21) and the resulting set of 
nonlinear coupled algebraic, equations was solved simultane- 
ously. The results did improve the prediction for the vertical 
heat flux ul and temperature variance f 2 but the comparison 
for radial heat flux and shear stress became worse. As noted 
by Gibson and Launder (1976), the above model is not good 
near an axis of symmetry. This is why, by incorporating them 
in the original model, no overall improvement in the prediction 
is achieved. 

Chen and Rodi (1975), Tamanini (1978), and Chen and Chen 
(1979) use the differential convective-transport equation (5) to 
determine z 2 in their predictions of buoyant jets. Equation (6) 


was used to eliminate e r Thus, the final form of the Z 2 equation 
becomes 



♦ .kS 

dr 


1 A 

r dr 



dr Rk 


(25) 


This equation was numerically solved for the temperature var- 
iance Z 2 with all other quantities needed to evaluate the coef- 
ficients determined from experiments. The value of C, was 
taken as 0.13. The best agreement, as shown in Fig. 7, was 
achieved with R = 0.35. When the average experimental value 
of R = 0.25 is used, the prediction peaks at about 6.0 (77 = 
0.04) as compared to the experimental Value of Z 2 of about 8.0 
(77 = 0.04). When the standard value of R = 0.8 is used, z 2 
is overpredicted by a factor of four. 

When Z 2 is calculated from the convective-transport equation 
(5), the diffusive transport is given by the simple gradient 
closure 

I (26) 

€ dr 


with C, = 0J 3 as given by Chen and Rodi (1980). The pre- 
diction for uz 2 , using Experimental information to evaluate the 
right-hand side of equation (26), is shown in Fig. 9. The pre- 
dicted curve peaks somewhat above and toward the centerline 
as compared to the data. 

Finally, we ask that if the models do not_depict the axial 
heat flux U7 and the temperature variance Z 2 correctly, then 
why do the predictions such as made by Hossain and Rodi 
(1982), Chen and Rodi (1975), Tamanini (1978), and Chen and 


Chen (1978) show reasonable agreement with the experiment 

for the mean velocity and buoyancy? The answer to this is 
that with R = 0.8 the temperature variance Z 2 from equation 
(21) or (25) is too large. This makes the vertical heat flux u 7 
from equation (19) large enough so that the mean velocity and 
buoyancy arc reasonably predicted. 

5 Summary and Conclusions 

The experimental data on buoyant plumes were used to 
evaluate various closure relations for turbulence transport. The 
objective was not to propose new models, but to evaluate the 
closure schemes proposed by other workers for buoyancy- 
dominated flows. The closures evaluated were those used in 
the k-c and algebraic stress models. The results are summarized 
below. 

1 The closure relations of the k-t model compare reason- 
ably with experimental data, except for the axial turbulent 
transport, which is drastically underprcdictcd. The axial heat 
flux governs the production due to buoyancy in the kinetic 
energy and dissipation equations and its correct prediction is 
very important. This is a probable reason why the results of 
Hossain and Rodi (1982) from the k-< model under predict the 
spreading rate for the plume by 10 percent even when axisym- 
metric jet corrections arc included. 

2 The ratio R of the time scales, which is used to determine 
the dissipation of the mean squared temperature in the alge- 
braic stress model, was found to be considerably different from 
the accepted value of R = 0.8. Apparently R is not a universal 
constant, but can vary from flow to flow and is influenced by 
the strength of the buoyancy present. From the experimental 
data on a plume it appears that R = 0.25 for strongly buoyant 
flows. 

3 The closure equations for the shear stress and radial heat 
flux of the algebraic stress models also compared well with 
experiment but are not better than the simple gradient closures 
used in the k-t model. The axial heat flux and mean squared 
temperature are predicted poorly in the central core of the flow 
and had incorrect trends. This drawback could be attributed 
to the assumption of local equilibrium, which resulted in the 
neglect of convection and diffusion terms in the transport 
equations for Reynolds stress and heat flux. However, no 
substantial improvement was. achieved by keeping the second- 
ary derivatives or by incorporating the model for the convec- 
tion and diffusion terms. Therefore, the full dynamic equations 
for Reynolds stress and heat flux with convection and diffusion 
arc required to predict the axial heat flux and temperature 
variance properly. 
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Abstract 

Equations for the mean and the turbulence quantities of compressible turbulent 
flows axe derived in this report. Both the conventional Reynolds average and the 
mass-weighted Favre average were employed to decompose the flow variable into a 
mean and a turbulent quantities. These equations axe to be used later in developing 
second-order Reynolds stress models for high-speed compressible flows. A few recent 
advances in modeling some of the terms in the equation due to compressibility effects 
axe also summarized. 
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Abstract 

Linear stability of flow in rotating pipe is studied. 
These flows depend on two parameters, which can be 
taken as the axial Reynolds number Re and the rotat- 
ing rate, q. In the region of Re 1 and q = 0(1), the 

most unstable modes are concentrated near the pipe 
wall, the so-called “wall” modes. These wall modes 
are found to satisfy a simpler set of equations contain- 
ing two parameters rather than four parameters as in 
the full linear stability problem. The set of equations 
is solved numerically and asymptotically over a wide 
range of the parameters. In the limit of Re — ♦ oo, 
the eigenvalue goes to the inviscid limit. The eigen- 
function shows a two layer structure. It reaches the 
inviscid limit over the main part of the domain, while 
near the wall of the pipe, the eigenfunction is repre- 
sented by a viscous solution of boundary layer type. 

1. Introduction 

Swirling flow is common in nature and technology. 
Fully-developed flow in a rotating pipe is an exact 
solution of the Navier-Stokes equations, and is the 
simplest available model of swirling flows. Swirling 
flows are known to be subject to instability, and the 
question of stability of flow in rotating pipe has con- 
sequently attracted a reasonable amount of attention. 
Pedley 1 showed that these flows are unstable to in- 
viscid non-axisymmetric perturbations when the ro- 
tation is fast (in a sense which will made definite). 
Stability to inviscid perturbations in the finite rota- 
tion rate region was studied numerically by Maslowe 2 . 
Later, Maslowe and Stewartson 3 extended this work 
in a significant way and established by asymptotic 
methods that the dominant unstable modes are wall 
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modes (that is, modes of motion concentrated asymp- 
totically close to the wall) in the limit of large az- 
imuthal wave number. 

For viscous perturbations, Pedley 4 , and simulta- 
neously Jesoph and Carmi 5 found that the critical 
Reynolds number at which the perturbation is neu- 
tral in the limit of fast rotation. Cotton and Salwen 6 
carried out comprehensive computations in search of 
neutral stability curves and they discovered that the 
neutral modes are center modes when Reynolds num- 
ber is large. Center modes in rotating pipe flow were 
later analyzed asymptotically by Stewartson, Ng and 
Brown 7 , and these authors speculated that center 
modes dominate for large Reynolds number. 

In this study, we investigate the effect of viscos- 
ity on the inviscid wall modes found by Maslowe and 
Stewartson 3 when the Reynolds number is large but 
finite. We find the proper scaling for Reynolds num- 
ber in order for viscous wall mode to exist and derive 
simplified governing equations for them. These equa- 
tions contain only two parameters instead of four pa- 
rameters in the full linear stability problem. The wall 
mode equations are then solved both numerically and 
asymptotically. 

The plan of this study is as follow: Linear stabil- 
ity analysis is formulated in section 2, where it is 
shown numerically that the most unstable modes are 
wall modes. The governing equations for viscous wadi 
modes are derived in section 3. Numerical solutions 
of the viscous wall mode equations are presented in 
section 4. An asymptotic analysis for the viscous wall 
modes equations is carried out in section 5, and sec- 
tion 6 concludes the paper. 

2. Linear Stability Formulation 

If length is nondimensionalized by the radius of the 
pipe L and velocity by the the axial velocity at the 
axis U , the laminar base flow in a pipe rotating with 
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angular velocity Cl is then described in a cylindrical 
(x,0, r) coordinate system by, 

U = (l — r’.gr.O) (1) 


where the inverse Rossby number 


9 = 


CIL 

IT 


( 2 ) 


measures the relative strength of rotation. This 
nondimensionalizaiion also defines a Reynolds num- 
ber 


Re = 


UL 


(3) 


where v is the kinematic viscosity of the fluid. 

Linear stability analysis concerns the stability of 
the motion subject to infinitesimal perturbations. 
Since the linear stability problem with base flow given 
by (1) is then separable in the x and 9 directions, the 
perturbation field may be written in the normal mode 
form. i.e. 


where 


7 = Jk(l — r 2 ) + m q — u 

The above equations for the perturbations are sup- 
plemented by the following boundary conditions. On 
the wall of the pipe, the no slip boundary conditions 
require 


u(l) = „(1) = tu(l) = 0 (5) 

At the center of the pipe, the perturbation must 
satisfy the following conditions to ensure that to be 
single- valued 8 

for m = 0 

u'(0) = v(0) = u>(0) = p'(0) = 0, 
for |m| = 1 

u(0) = i/(0) + ini u/(0) = p(0) = 0, 
for m otherwise 

u(0) = i/(0) = n/(0) = p( 0) = 0. (6) 


A [u(r), v(r), w(t), p{r)]exp [i(fcx + m 3 - art)] 


where A is an arbitrary constant, k is the axial wave 
number of the perturbation and m is the azimuthal 
wave number of the perturbation. u> is the complex 
frequency, with its real part being the frequency and 
the imaginary part being the growth rate. Without 
loss of generality, m is taken as positive, k is any real 
number, and u/ is to be found. If Im(o/) is positive, 
the flow is linearly unstable. 

The Navier-Stokes equations linearized about the 
base flow, and the equation of continuity are of the 
following form in the cylindrical coordinate system 
used, 


1*711 — 2 rw + xkp = 

1 / d 2 u 1 du 

Re \dr 2 + r dr 



imp 

171/ + 2 qw -1 = 

r 

1 / d 2 v 1 dv 

Re V dr 2 ^ r dr 


m 2 + 1 
+ 


2 177111/ 



dp 

ijw - 2 qv + — = 


1 / d 2 w 1 dw 
Re \ dr 2 ^ r dr 


771 2 + 1 

1 — u/ 


2 imv 

r 2 



•T l7Tl 

iku + v 

r 


dw w 



= 0 


(4) 


The above ordinary differential equations (4) and 
the boundary conditions (5) and (6) form an eigen- 
value problem with w as the eigenvalue. Nontrivial 
solution exists only when u takes some specific values 
given by 


w = v{Re % q; m,k) 

Iflm(u/) is less than zero, the flow is stable; if Im(u/) 
is positive, the flow is said to be linearly unstable. We 
are interested in the regions where the instability oc- 
curs. Earlier studies show, and our results confirm, 
that instability occurs only for mk < 0. In the fol- 
lowing, we will take m positive, and k negative. 

The linear stability problem was studied exten- 
sively by numerical means by Cotton and Salwen® 
and by Yang®. Cotton and Salwen searched in the 
parameter space (Re, q ; m, Jk) for the neutral modes 
(modes with zero o/ t ), and found that the neutral 
modes are center modes, with most of the nontrivial 
activities confined near the center of the pipe. Yang’s 
study concentrated on the most unstable modes, and 
we shall briefly describe the most unstable modes that 
he found in the region Re 1 and q = 0(1). 

Table 1 gives the information for the most unsta- 
ble modes at Re = 10000 and some different values 
of q 1 s. As q is increased, the the most unstable mode 
has larger values of both the azimuthal and axial wave 
numbers. In addition, we find for q = 0(1), although 
the real part and the imaginary part of the eigen- 
value have quite different size, the difference of the 
real part of the eigenvalue from mq is of the same 
order as that of the imaginary part. These findings 
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Tabic 1: The most unstable modes for different q's 
with Re = 10000. 

m: the azimuthal wave number, 
k : the axial wave number, 

(j: the eigenvalue of the most unstable mode. 


<? 

771 

k 

u> 

0.5 

3 

-0.50 

(1.25, 0.177) 

1.0 

5 

-0.81 

(4.67, 0.296) 

1.5 

7 

-1.04 

(10.13, 0.379) 

2.0 

9 

-1.20 

(17.61, 0.439) 

3.0 

11 

-1.22 

(32.63, 0.518) 


suggest the viscous wall mode scalings studied in the 
next section. For q = 3, the most unstable mode has 
the wave numbers m = 11, Jk = —1.22. The eigen- 
function corresponding to this most unstable mode is 
shown in Fig 1. The eigenfunction is normalized such 
that the maximum axial velocity is 1. The real part 
of the eigenfunction is drawn in solid lines, and the 
imaginary part is drawn in dotted lines. The non- 
trivial behavior of this eigenfunction takes place in a 
thin region near the wall of the pipe. This behavior 
justifies the “wall mode” terminology used, and will 
be the subject of further study in next section. 


u 


M 



Figure 1: Eigenfunction for Re = 10000, q = 3, m = 
11, & = —1.22. (a) axial velocity, (b) azimuthal ve- 
locity, (c) radial velocity. 


3. The viscous wall mode equation 

From the numerical computations in the last sec- 
tion, it is clear that in the region of Re 1 and 
q = 0(1), the dominant modes are given by the 
asymptotic wall modes, a type of modes with large 
azimuthal wave number and with nontrivial behavior 
concentrated near the wall. 

When azimuthal wave number m is large, there 
could be another type of mode for general swirling 
flows, the ring mode, as demonstrated by Leibovich 
and Stewartson 10 . In the case of rotating pipe flow, 
ring modes do not exist. Because of the existence of 
the pipe wall, wall modes characterize the behavior 
of perturbation with large azimuthal wave number. 

Maslowe and Stewartson 3 studied the stability of 
inviscid pipe flow to perturbations of very large az- 
imuthal wave number, and established that the pre- 
vailing modes are the wall modes in the limit of 
m — ► oo. Stewartson 11 analyzed the effect of the 
viscosity on ring modes and found that viscous ef- 
fects come into play in higher order terms because 
the inviscid ring mode solution can satisfy the exact 
boundary conditions of the problem. 

We study the wall mode when the viscosity is taken 
into account. Since the inviscid wall mode is close to 
the wall, where the no-slip condition is required but 


is not satisfied by the inviscid solution, it is expected 
that the viscous effect might be important in this 
case, at least in certain regions . 

We are going to study the linear stability problem 
for large azimuthal wave number and large Reynolds 
number, i.e. both m and Re are large. In this case, 
the proper wall mode scalings are 


(m 3 + k 2 ) 1 / 2 

w = iw 

? = % 

771 

r . , = -fetel ±££!1 

2k 
2k Q 

H (m 3 + Jk 3 ) 1 / 3 

The scaling for the radial variable means that the be- 
havior of the wall mode is confined to a small distance 
to the wall comparable with the wave length of the 
perturbation. The scaling for the Reynolds number 
giv^s the balance between the viscous term and the 
inertia term. The form of the scaling for the complex 
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frequency is suggested by the numerical computations 
of the full linear stability problem. In above, a mi- 
nus sign is introduced in places where k appears, for 
we know from the full linear stability problem that 
instability occurs only when mk < 0. 

Upon substituting the above expressions into the 
linearized Navier-Stokes equations and the equation 
of continuity for the perturbation, and dropping the 
terms of order m“ l and smaller, we find that w should 
satisfy the following single equation, after u, v y p are 
eliminated. 


nonlinear eigenvalue problem. This nonlinear eigen- 
value problem is changed to a system of linear eigen- 
value equations by letting 

Y = Lw (12) 

In Y } w , the governing equations are 

L(D 2 -l)Y + 2 DY -(Q - 2 )w = 0 

Y -Lw = 0 ( 13 ) 


LDLDw — LLw -f LDw — Qw ~ 0 (8) 


where 


D = — 


d 

drj 


and 


L = (£)’ -!)-(„ + *) 

Re 


g(m + qfc) 
W k 


( 9 ) 


which emerges as one of the independent parameters 
for the viscous wall mode. 

The boundary conditions are 


w = Dw = D 2 Lw — Lw = 0 at 77 = 0 (10) 

and 

w = Dw = D 2 Lw — Lw = 0 as 77 — ► 00 ( 11 ) 

Thus, we have established the formulation for the 
viscous wall modes. The governing equation takes 
a simpler form compared with the full problem, and 
the number of the independent parameters is reduced 
from four to two. Of these two parameters, Re and Q } 
Q measures the effect of rotation on the wall mode 
and Re measures the effect of viscosity on the wall 
mode. The solutions of above wall mode problem 
will give: 

w = w{Re, Q) 

If Im(ui) is less than zero, the flow is linearly stable; 
if Im(cD) is positive, the flow is said to be linearly 
unstable. 


4 . Numerical Solution of Wall Mode Equations 

In general, solutions of the wall mode equation have 
to be found numerically. Because the eigenvalue en- 
ters quadratically, this wall mode equation gives a 


and the boundary conditions are 

w = Dw = D 2 Y — Y = 0 at 77 = 0 (14) 

w = Dw = D 2 Y — Y = 0 as 77 00 ( 15 ) 

We employed a spectral method with Chebyshev 
polynomials as the basis functions to solve equations 
( 13 ) - ( 15 ). Because the domain of definition extends 
to infinity, while Chebyshev polynomials are only de- 
fined over [-1, l], the method of domain truncation 
was used to numerically truncate the domain of def- 
inition from [ 0 , 00) to [ 0 , b] and the boundary condi- 
tions at infinity are replaced by 

w — Dw = D 2 Y — Y = 0 at 77 = 6 ( 16 ) 

In the spectral method used, we write 
N 

w = ^w,Ti_ l (y) 

»=1 

AT 

Y = ( 17 ) 

»=1 

where y is related to 77 by 

y = 277/6 — 1 . 

Thus the domain of definition for y is [-1,1]. The 
reduction from the differential equations to a set of 
algebraic equations is made by a Galerkin-Tau projec- 
tion — i.e. the Galerkin method is used to project the 
equations while the Tau method is used to enforce the 
boundary conditions on the spectral representations 
of the perturbation field. The Galerldn-Tau projec- 
tion results in an set of algebraic equations, which are 
of the form of generalized eigenvalue problem with 
complex matrices. 

Two methods are used in this study to find the 
eigenvalues and the eigenvectors of this complex gen- 
eralized eigenvalue problem. One method uses the 
IMSL subroutine EIGZC which uses the QZ transfor- 
mation to find all the eigenvalues and, optionally, all 
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the eigenvectors. The other method used is an inverse 
power iteration for the generalized eigenvalue prob- 
lem developed by Kribus 12 , which finds the eigen- 
value closest to the initial guess and its corresponding 
eigenvector. The inverse power iteration is faster, so 
it is used whenever a good guess is available. The QZ 
is used to provide the starting values for the Inverse 
Power Iteration. It is also used when the phenomenon 
of mode jumping is suspected to occur. 

Since the eigenfunctions decay exponentially as 
77 — * oo, b = 10 was found sufficient for most of the 
calculations. The number of terms in the Chebyshev 
representation, N , varies with Re. For order one Re , 
we must take N = 60 , and b = 10 to achieve three 
digit accuracy in eigenvalue, and O(10 -4 ) accuracy 
for the eigenfunction. But when Re is large, the so- 
lution shows a behavior of boundary type for 77 near 
zero, i.e. near the wall of the pipe, and a large value 
of AT is needed to resolved this region. The largest N 
needed for the parameter range covered here is 115 . 

There is a symmetry in the eigenvalue problem due 
to the replacement of the boundary conditions at in- 
finity by the conditions at 77 = b. The solution is 
invariant under 

77 «— ► 6 — 77 

Q «— ► Q* — b 

w «— ► w* 

Y ~ Y m 

where the star means the complex conjugate. This 
symmetry signifies that for a given solution, its im- 
age about 77 = 6/2 is also a solution of this equation 
with the same growth rate. Apparently, this solution 
is a spurious mode in the sense that it is a solution 
of the differential equation after enforcing the bound- 
ary condition at finite b rather than the solution of 
the original differential equation defined over the in- 
finite domain. This symmetry property can be used 
to check the resolution of the numerical solution of 
the algebraic problem, which should also exhibit this 
symmetry. 

Extensive computations were carried out in the 
(Re, Q) plane. Fig 2 shows the imaginary part of the 
eigenvalue u>, which is proportional to the growth rate 
of the perturbation, vs. Re for some different values 
of Q. The eigenvalues shown are for Q = 5 , 10, 20 , 30 , 
40 , 50 respectively, although numerical computations 
were carried out for a larger range of Q. The real 
part of the eigenvalue is shown in Fig 3 . The growth 
rate increases as Q is increased, which ni^ans that ro- 
tation helps perturbations to extract energy from the 



Figure 2: The imaginary part of the wall mode eigen- 
values vs. Re for Q = 5 , 10 , 20 , 30 , 40 , 50 . 


base flow. This result is in agreement with that of 
Pedley 1 , who found that the maximum growth rate 
is 2.0 and is reached in the fast rotation limi t. This is 
exactly the same as the upper bound for the growth 
rate for flow in the rotating pipe, as shown by Joseph 
and Carmi 5 . 

The growth rate increases with Re. As Re gets 
large, the growth rate will increase. The eigenvalues 
seem to change smoothly as Re — ► 00 and approach 
to the results found by Maslowe and Stewartson from 
their inviscid analysis. The limit of Re — ♦ 00 will be 
analyzed asymptotically in the next section, and a 
comparison of the results with the inviscid case will 
be made. 

Fig 4 shows the eigenfunction for Re = 2 , and Q = 
5 . The eigenfunction plotted is normalized such that 
its maximum modulus is 1 and its phase at the posi- 
tion of maximum modulus is 0 . As in the eigenfunc- 
tion plotted in Fig 1, the real part of the eigenfunction 
is drawn in solid line while the imaginary part of the 
eigenvalue is drawn in dotted line. To see the effect 
of increasing Q on the eigenfunction, we show in Fig 
5 the eigenfunction for Re = 2 and Q = 20 . As Q 
increases, the eigenfunction is pushed outward, but 
this response is not very sensitive to Q. To see the 
effect of changing Re, Fig 6 shows the eigenfunction 
for Re = 100 and Q = 5 . As Re increases, the eigen- 
functions are pushed towards 77 = 0, i.e. toward the 
wall. The limit of Re — * 00 will be the further studied 
in the next section. 

5 . The Limit of Re — » 00 
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Figure 3: The real part of the wall mode eigenvalues Figure 5: Wall mode eigenfunction, Re = 2, Q = 20 
vs. Re for Q = 5, 10, 20, 30, 40, 50. 




Figure 4: Wall mode eigenfunction, Re = 2, Q = 5 Figure 6: Wall mode eigenfunction, Re = 100, Q = 5. 
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Maslowe and Stewartson 3 carried out the wall 
mode analysis for m 1 for the inviscid case. One 
of the purposes of our viscous wall mode analysis is 
to see if the results of their inviscid analysis is the 
limit of the viscous analysis for large Reynolds num- 
ber. We use a perturbation technique to address this 
question in this section. 

For Re ^ 1, w can be expanded formally by taking 
Re 1 as the small parameter, we write: 

u> = WS{n) + Re~ l W?(T,) + ... 

To leading order, the equation for Wg is: 

° lw o - U + = 0 ( 18 ) 

(77 + Wj' 

with the following boundary conditions: 

W 0 °( 0) = 0 (19) 

WS( 77-00) - 0 (20) 

This poses the same eigenvalue problem that was 
studied by Maslowe and Stewartson, thus their solu- 
tions (both the eigenvalue and the eigenfunction) may 
be viewed, as might have been expected, as the first 
term in a formal outer expansion in inverse powers of 
Re. 

The outer solution thus found can satisfy all the 
boundary conditions at infinity since the outer solu- 
tion decays exponentially as 77 — ♦ 00. But not all the 
boundary conditions at 77 = 0 can be satisfied. For 
example, DWg( 0) 0. Thus, another (inner) solu- 

tion of boundary type near 77 = 0 is needed. 

The inner variable is found to be: 

< = nS*'' 2 

The inner expansion is assumed to be: 

* = Re d (Wo«) + Re~ l,7 Wi( 0 + ...) 

where d is to be determined. 

After substituted into the governing equation, the 
equation to leading order is found to be: 

(D 2 +u) 2 D 2 Wi = 0 (21) 

where 



The boundary conditions for the inner solution are: 

0 ) = 0 

DW*( 0) = 0 

(D 1 + u,) 2 W<( 0 ) = 0 (22) 


In addition, the solution is required to match the 
outer solution in the matched asymptotic sense. 

The general solutions of the inner problem in the 
leading order can be found, and they are 

wt - Pi exp(AiC) + Pi < exp(AjC) + P 3 exp(AjC) 
+0* C er p(-^jC) + Pi + P« C (23) 

where 

Ai ,2 = 

and pi, Pi, Pz, P4, Ps, /?© are constants to be deter- 
mined by the boundary conditions and the matching 
conditions, u is given by the outer solution. 

The boundary conditions at £ = 0 require that 

Pi + Pz 4- Ps = 0 

A^i 4- Pi 4* PzXi 4- P\ 4- Pe = 0 

2A J Pi 4- 2A 3 P\ = 0 

As £ — ► oo, the inner solution must also match the 
outer solution. Near 77 = 0, the outer solution is, to 
leading order, 

Wg = a 77 

where a = DWg( 0) 9^ 0. 

To carry out the matching, we need to know the 
behavior of Wq for large £. The behaviors of the 
exponential terms are determined by the sign of the 
real parts of Ai, A 2. Since Ai, A? are the square roots 
of (ia>), Ax is the negative of A 3 . In this study, we 
take 

Re{ A x ) < 0 

Re{Xi) > 0 

As C — ♦ 00, the matching of the exponentially grow- 
ing terms gives: 

Pz = 0 
Pa = 0 

The exponentially small terms are immaterial, and 
the matching of the algebraic terms gives 

Re d Pe( = <xv (24) 

which gives 
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and 


0e = cr 

All the other constants can be determined, yielding 
0 \ = — <*/Ai 

02 = 0 
0s — ct/Ai 

The inner solution, to the leading order, is 

WS(0 = o(C + i-(l-«p(A 1 0)) (25) 

*1 

Thus, we have the following composite solution to the 
leading order: 

WS = wg( v ) + W5(0 - av (26) 

Thus, one sees that indeed, the solutions for the 
most unstable modes found by Maslowe and Stew- 
artson are the correct limit when Re — ► oo, except in 
a thin layer near the wall of the pipe where the flow 
field is described by a viscous layer of the boundary 
layer type. The eigenvalues found are the same as for 
the inviscid case. 

The above asymptotic analysis is confirmed by the 
numerical calculation. In Table 2, we present the 
eigenvalues for Re = 100 and in Table 3, we present 
the eigenvalues from the inviscid calculation for,a few 
values of Q. It is seen that differences of the eigen- 
values in these two cases are small. In Fig 7, we show 
the eigenfunction for Q = 5.0 from the inviscid calcu- 
lation, It is seen that the general shape agrees with 
the viscous calculation presented in Fig 6. To see 
the existence of a thin viscous layer near the wall, we 
present in Fig 8a a blow-up of Fig 6, and in Fig 8b 
a blow up of Fig 7 near the wall. It is seen that flow 
fields near the wall are different, the viscous solution 
has a zero slope while the solution from the inviscid 
equation has a non-zero slope. 

6. Discussion 

We have examined the linear stability of rotating 
pipe flow to perturbations of large azimuthal wave 
number. It is found that when the azimuthal wave 
number is large, the nontrivial behaviors are concen- 
trated near the wall of the pipe, so the prevailing 
variations are manifested as wall modes. The equa- 
tions governing wall modes are found to contain two 
parameters Re and Q, which measure the Reynolds 
number and swirl, respectively. 


Table 2: Eigenvalue of viscous wall mode. 


Q 

eigenvalue (Re = 100) 

5.0 

(-0.14917D+01, 0.72752D+00) 

10.0 

(-0.18163D+01, 0.15096D+01) 

20.0 

(-0.21617D+01, 0.26668D+01) 

30.0 

(-0.23766D+01, 0.35755D+01) 

40.0 

(-0.25358D+01, 0.43505D+01) 

50.0 

(-0.26633D+01, 0.50384D+01) 


Table 3: Eigenvalue of inviscid wall mode. 


Q 

eigenvalue (Re = oo) 

5.0 

(-0.14142D+01, 0.68305D+00) 

10.0 

(-0.17476D+01, 0.14875D+01) 

20.0 

(-0.21038D+01, 0.26574D+01) 

30.0 

(-0.23247D+01, 0.35710D+01) 

40.0 

(-0.24877D+01, 0.43487D+01) 

50.0 

(-0.26181D+01, 0.50383D+01) 



Figure 7: Wail mode eigenfunction for the inviscid 
case, Q = 5. 
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Figure 8: Comparison of the wall mode eigenfunc- 
tions near the wall, (a) the viscous solution at Re = 
100, (b) the inviscid solution. 

For large azimuthal wave number, the wall mode 
is a distance 0(m -1 ) away from the wall. Vis- 
cous effects are confined to a layer with thickness 
of 0(Re ” 1 / 2 ) near the wall. The case treated here 
is when those two layers are of the same order, i.e. 
Re = 0(m 2 ). When Re/m 2 1, we would ex- 

pect the solution to be mainly inviscid, as found by 
Maslowe and Stewartson, except in a thin viscous 
layer near the wall. This is indeed the case as shown 
both numerically and asymptotically. 

The region for the rotation rate considered is q ■= 
0(1). For large value of q , which corresponds to the 
fast rotation case, Pedley 4 shows that as Reynolds 
number is increased, the dominant mode is taken by 
perturbations of larger and larger values of m, the 
azimuthal wave number. In the limit of Re — ♦ oo, the 
azimuthal wave number for the dominant modes will 
also go to infinity, and these modes could be viewed 
as wall modes. Thus, our work here can be viewed 
as an extension of Pedley’s work 1,4 to the region of 
, = 0 ( 1 ). 

Left unexamined is the case when the inviscid so- 
lution has a singularity, in which case the form of 
the expansion we assumed would break down. In the 
study of Maslowe and Stewartson, the leading mode 
is a neutral mode when Q = 2. For Q > 2, the 
leading mode is unstable, but there are still neutral 
modes. The neutral modes have a l/(r — r 0 ) singu- 
larity, where r 0 is the location of the singularity. It is 


expected that a local viscous critical layer based on 
our viscous wall mode formulation would get rid of 
this singularity. 
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Development of A Nev,’ Flux Splitting Scheme 


Meng-Sing Liou* and Christopher J. Steffen, Jr.t 
NASA Lewis Research Center 
Cleveland, Ohio 


Maximizing both accuracy and efficiency has been the 
primary objective in designing a numerical algorithm 
for computational fluid dynamics (CFD). This is espe- 
cially important for solution of complex 3D problems 
which often involve Navier-Stokes equations with turbu- 
lence modeling and chemical species equations. Upwind 
schemes have been well received for both their capa- 
bility of resolving discontinuities and their sound theo- 
retical basis in characteristic theory for hyperbolic sys- 
tems. Several flux splitting schemes, notably by Steger- 
Warming, Van Leer, Osher, and Roe, have been tested 
and discussed extensively in the past decade. 

However, several inherent shortcomings exist in each 
of these schemes. For example, while the Van Leer 
scheme is simple, taking only 0(n ) operations, n being 
the number of equations, and yields accurate solutions 
for inviscid problems, it suffers inaccuracy for predicting 
velocity and temperature fields in the viscous problems. 
The Roe scheme, commonly accepted as the most accu- 
rate scheme available currently, however is a great deal 
more complex and costly due to the matrix operation, 
requiring 0(n 2 ) operations. Moreover, the extension to 
the chemically reacting flows renders no unique way of 
defining the ‘Roe-averaged* states. The Osher scheme 
has the smooth property and has recently been general- 
ized by Suresh and Liou to deal with chemically reacting 
flows. But the determination of the intermediate states 
also requires 0(n 2 ) operations. Thus it is logical to ask 
whether there is room in the universe of upwind schemes 
for improvement to arrive at a simple (O(n) operations) 
and accurate scheme for a wide range of problems. 

In this paper, we summarize recent successes of a 
new splitting scheme for some model aerodynamic prob- 
lems where Van Leer and Roe schemes failed. The new 
scheme is based on a rather different idea of splitting in 
which the convective and pressure terms are separated 
and treated differently in accordance with the under- 
lying physical intuitions. We propose an appropriately 
defined cell-face advection Mach number using values 
from the two strawidling cells via associated characteris- 
tic speeds. This interface Mach number is then used to 
determine the upwind extrapolation for the convective 
quantities. Next the pressure splitting is weighted us- 
ing polynomial expansions of the characteristic speeds. 


Thus, the name of the present scheme is properly coined 
as Advection Upstream Splitting Method (AUSM). The 
scheme is remarkably simple and yet its accuracy in the 
present study rivals and in some cases surpasses the Roe 
scheme in the Euler and Navier-Stokes solutions at con- 
siderably reduced computational effort. The detailed 
formulation of the scheme is net shown here due to space 
limitation. However it will appear elsewhere. 

The calculation of the hypersonic conical flow demon- 
strates the accuracy of the splittings in resolving the 
flow in the presence of strong gradients. The tempera- 
ture and pressure profiles of the first order results are 
shown in Figs. 1 (a) and (b). The Van Leer splitting is 
seen to produce a thicker boundary layer which in turn 
further displaces the shock wave. Both AUSM and Roe 
solutions are in excellent agreement. 

The second series of tests involve the 2D inviscid flow 
over a NACA 0012 airfoil. The results, not included 
here, demonstrate that the level of entropy generation 
at the stagnation point is about three times smaller than 
the Roe solution. 

In the third case we calculate a series of supersonic 
flows over a circular cylinder. The Roe splitting in all 
conditions and grids tested yields anomalous solutions 
(sometimes referred to as the carbuncle phenomenon), 
which may appear as nonsymmetric, protuberant, or in- 
dented contours, see Fig. 2. The mode of these non- 
physical solutions appears to be sensitive to changes in 
Mach number or grid. The AUSM, however, gives ex- 
pected solutions in all calculations. 

The fourth test deals with a 2D shock wave/laminar 
boundary-layer interaction. In Fig. 3, the AUSM is seen 
to give excellent agreement with the data especially in 
the reattachment region, which has been in defiance of 
many previous calculations in the literature. Also the 
oblique shock appears to be more tightly captured by 
the AUSM. 

In summary, it appears that the new splitting scheme, 
AUSM, has delivered the promise by improving the ac- 
curacy as well as efficiency significantly. As usual, the 
final judgment will be decided via many more tests and 
further modifications of the scheme. 
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Fig. 1 Conic flow with — 7.95, Re = 4.2 x 10 5 ,and half cone angle = 10°: (a) pressure, and (b) temperature 



Fig. 2 Mach contours for a Mach 6 inviscid flow over a circular cylinder: (a) the AUSM solution, (b) the Roe 
solution displaying a protuberrant, two-shock contours. 



Fig. 3 Mach contours and comparison of skin friction with data for M Q Q = 2.0, Re = 2.96 x 10 5 [Hakkinen et al]: 
(a) the AUSM solution, and (b) the Roe solution. 145 


High-Order Polynomial Expansions(HOPE) for Flux-Vector Splitting 


Meng-Sing Liouf and Chris J. Steffen, Jr.* * 


NASA Lewis Research Center 
Cleveland, OH 44135, U.S.A. 


Summary 

The Van Leer flux splitting is known to produce excessive numerical dissipation for 
Navier-Stokes calculations. One example is the incorrect prediction of boundary-layer 
profiles. We attempt in this paper to remedy this deficiency by introducing a higher- 
order polynomial expansion(HOPE for short) for the mass flux. In addition to Van 
Leer’s splitting, a term is introduced so that the mass diffusion error vanishes at M = 0. 
Several splittings for pressure are proposed and examined. The effectiveness of the 
HOPE scheme is illustrated for 1-D hypersonic conical viscous flow and 2-D supersonic 
shock-wave/boundary-layer interactions. Also, we give the weakness and suggest areas 
for further investigation of the scheme. 


Introduction 


In the past decade, upwind differencing schemes have gained considerable attention 
for their accuracy and robustness in Euler flows with discontinuities, shock waves in 
particular. Naturally, significant research effort in the CFD community has been focused 
on maxmizing the accuracy and efficiency, among other objectives. Four popular but 
conceptually different flux splitting ideas have been utilized for nearly 10 years: Steger 
and Warming, Van Leer, Roe, and Osher. However, each scheme has an associated 
weakness when numerical accuracy and efficiency are considered. 

In this paper, we deal specifically with the improvement of Van Leer’s flux vector split- 
ting[l]. Besides its simplicity, Van Leer’s splitting has the following properties: (1) it can 
be interpreted as a special member of a family of second-order polynomial expansions[2], 
and (2) the associated flux Jacobian and eigenvalues are continuous at the sonic points. 
Van Leer’s choice allows one vanishing eigenvalue in the case of an ideal gas, thereby 
resulting in a crisp shock representation. Furthermore, the continuous differentiability 
is helpful for convergence acceleration, e.g., in multigrid schemes. 

However, failing to recognize the contact discontinuity, the Vain Leer splitting[l] pro- 
duces excessive numerical diffusion and thus requires a huge number of cells to correctly 
resolve the boundary-layer flow. Some improvements have been demonstrated recently 
by Hanel et al[2] and Van Leer[3] for 1-D conical, hypersonic viscous flow, but a pressure 
glitch arises. A new scheme by the present authors[4] has been proposed that not only 
corrects this pressure difficulty, but also is remarkably simple to implement. Neverthe- 
less, the above schemes[2-4] have already departed from the ideas of flux vector splitting 
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and in fact become more like the flux difference splitting. Since the differentiability and 
simplicity are desirable properties, one would still wish to search for a better splitting 
scheme that is strictly based on the flux vector splitting. 

In this paper, we propose a family of higher-order polynomial expansions for the mass 
flux that diminishes the diffusion error as M — ► 0. We give a detailed study of the 
accuracy of the scheme for 1-D conical flow and 2-D shock wave/boundary-layer inter- 
actions. The weakness of the scheme is also pointed out and possibile improvements 
suggested. 


Analysis 

To exemplify the concept, let us consider the quasi two-dimensional system of equations 
for conical flows: 

8V 8F 
dt + d r? “ s 

where U T = (p, ptx, pv, pE), F T = (pv,pvu, pv 2 + p,pvH ), E = e + 1/2 (u 2 + v 2 ), and 
H = E + p/p. The flow considered consists of a very thin shear layer at the wall 
and a shock wave away from the wall. An algorithm must be capable of minimizing 
the numerical smearing(diffusion) at the locations where an eigenvalue changes sign or 
approaches zero. For example, Van Leer’s splitting[l] can represent shock profile well, 
while greatly diffusing the boundary layer. The Van Leer split mass fluxes are: 

F i = F+ + F~\ Ft = ±pa/4(M ± l) 2 . 

The net difference from the curve it approximates is largest at Af = 0; its value equals 
pa/2. This error, viz numerical diffusion, significantly broadens the boundary layer, 
leading to incorrect velocity and temperature profiles. A simple way to remove the 
diffusion at Af = 0 is by adding an extra higher-order term that allows the split mass 
fluxes to pass through the origin(Fig. 1), i.e., 

Ft = ±pa/4[(M ± l) 2 + m 1 (M)(M 2 - l) 2 ], 

where the higher-order term has a coefficient mi, in general function of Af. It should 
have the following properties: 

(1) mi — ► — 1 as Af — ♦ 0; 

(2) m,{M) = 

(3) mi — ► 0 as Af — ► ±1. 

A formula satisfying those properties is chosen as: 

fm * (Af a - !)/(Af a + i) s , 

where the exponent 5 is a free parameter; also shown in Fig. 1 is mi vs M with 5 = 4. 
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Fig. 1. Mass Flux Splitting in HOPE. Fig. 2. Pressure Splitting. 


In the conical flow calculations, the accuracy and convergence appear to be insensitive 
to the specified values of 5 = 2,4,6. Now, regarding the flux as a sum of convective and 
pressure terms, we can write the splitting formula for the flux vector: 



/ pv \ ± 

pvu = F ± 
pvv + p \ 1 

' pvH > 


1 \ 

u 

t; 

' h' 


+ 


/ o \ 
0 



With the realization in [5] that the pressure splitting could be considered separately in 
the Van Leer formula[l], a whole host of freedoms for the pressure splitting becomes 
possible. Following is the list of formulas tested: 


(Pl): 

p ± = ?l/4(M±l) 3 (MT2)p, 

(P 2): 

P* = (pl)=F3/4M(M 3 -l) J p, 

(p3) •' 

p k = (pi) ± 3/4mi M (M 3 — l) 3 p, 

and 


(p4) : 

p* = 1/2(1 ±7M)p. 

Figure 2 displays the distribution of the split pressure vs M. The first formula is that 
used by Van Leer[l]. The second and third splits, (p2) and (p3), yield vanishing 


slope at M = 0, thus corresponding to central differencing. However, no instability 
was encountered in the conical flow problem with the (p2) or (p3) split used in an 
implicit code. The fourth split (p4) is obtained from an approximate integration along 
characterics. As will be seen later, the four formulas give essentially the same results 
for the conical flow calculated. 


Results And Discussion 


In this paper, two cases were tested to check the accuracy and convergence of the 
HOPE scheme. The first case is the 1-D self-similar conical flow over a 10-degree 
half cone at hypersonic speed, for which a detailed comparison study was conducted. 
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The flow conditions are: M 0 0 = 7.95, and Re = 4.2 x 10 5 . Since P r = 1.0, exact 
solution gives adiabatic wall temperature, 13.64TOQ. The second case is the 2-D shock 
wave/laminar boundary-layer interactions, for which experimental measurements were 
available[6j. The conditions are: M <» = 2.0, .Re*, = 2.96 X 10 5 , and oblique shock 
angle (3 = 32.585 degrees. In both cases, the results from the Roe splitting are also 
included for comparison. An implicit Newton iteration procedure was used to achieve 
steady-state solution with L 0 0 residual dropped by five orders of magnitude. 

Figures 3 and 4 show the pressure and temperature distributions from the first- and 
second-order solution on a 65-grid; little difference is seen. A monotone solution across 
the shock is obtained with the first-order scheme while oscillation appears in the second- 
order scheme, which can be eliminated by a TVD procedure. It is noted that the first- 
order pressure is smooth at the edge of the boundary layer, unlike the Roe solution which 
shows a slight discontinuity(not shown here). Although the boundary layer exhibits a 
steep temperature gradient, the HOPE scheme predicts the wall temperature correctly, 
indicating removal of the numerical diffusion associated with the original Van Leer 
splitting. 
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Fig. 3. Pressure Profile of Conic Flow. Fig. 4. Temperature Profile of Conic Flow. 


Figure 5 displays the results using various pressure splittings; they are practically iden- 
tical except the Van Leer pressure split (pi) shows some minor oscillation near the wall. 
However, the pressure splittings show significant effect on the convergence rate. The 
(p3) and (p4) splits are the best, comparable to the Roe splitting, while the other two 
are roughly two to three times slower. These may indicate possible instability in a more 
complex case. 

Finally, for the 2-D case, the surface pressure and friction coefficient are plotted in 
Figs. 7 and 8. The first-order HOPE results compare fairly with Roe’s splitting and 
experimental data. However, the second-order calculation experienced difficulty in con- 
vergence in which the residual was reduced by only two orders of magnitude and the 
result is not presented here. 
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Fig. 5. Comparison of Pressure Splits. Fig. 6. Convergence History. 


We suspect that a further investigation on other pressure splittings may lead to success 
in stability and convergence. Nevertheless, a systematic study of the eigenvalues of 
the split fluxes and the complete discretized system will prove to be a useful endeavor. 
Above all, the present research suggests that there are still possibilities in flux-vector 
splitting after Van Leer’s appeared nearly 10 years ago. The possibilities may very well 
still lie in the mass-flux splitting. 




Fig. 8. Friction Coefficient at the Wall. 
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Using an LU Scheme 
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ABSTRACT 

A new three-dimensional numerical program that incorporates comprehensive real 
gas property models has been developed to simulate supersonic reacting flows. The code 
employs an implicit, finite volume, Lower-Upper (LU), time-marching method to solve 
the complete Navier-Stokes and species equations in a fully-coupled and very efficient 
manner. A chemistry model with nine species and eighteen reaction steps is adopted 
in the program to represent the chemical reactions of H2 and air. To demonstrate the 
capability of the program, flow fields of underexpanded hydrogen jets transversely injected 
into the supersonic airstream inside the combustors of scramjets are calculated. Results 
clearly depict the flow characteristics, including the shock structure, the separated flow 
regions around the injector, and the distribution of the combustion products. 
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ABSTRACT 

A two-dimensional, compressible Navier-Stokes 
equations with a k — c turbulence model are solved 
numerically to simulate the flows of compressible free 
shear layers. The appropriate form of k and e equations 
for compressible flows are discussed. Sarkar’s modelling 
is adopted to simulate the compressibility effects in the 
k and e equations. The numerical results show that the 
the spreading rate of the shear layers decreases with 
increasing convective Mach number. In addition, fa- 
vorable comparison was found between the calculated 
results and Goebel and Dutton’s experimental data. 


INTRODUCTION 

Recent national interest in trans-atmospheric ve- 
hicle has rekindled the hypersonic research. For this 
vehicle, a supersonic combustion ramjet (scramjet) en- 
gine was proposed to provide the power. Inside this 
scramjet engine, compressible mixing layers axe impor- 
tant phenomena. The performance of the engine will 
depend on the supersonic mixing and the flame holding 
of shear layers. 

The behavior of incompressible mixing layers has 
been studied extensively. However, additional study 
is required to understand the compressibility effects of 
free shear layers at high speeds. Figure 1 shows a sketch 
of a typical free shear layer. Two streams at different 
temperatures, densities, and Mach numbers merge to- 
gether to form a free shear layer. Various combinations 
of flow conditions of high-speed streams and low-speed 
streams allow for the systematic study 
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of compressible shear layers. In this paper, we report 
the incorporation of a k — e model with compressibility 
effects to a two-dimensional flow equations solver for 
the simulation of compressible shear layers. First, we 
point out the derivation procedure of the compressible 
k and c equations. Unlike the procedure in the incom- 
pressible equations, both Favre and Reynolds averaging 
procedures 1 are performed in deriving the flow and tur- 
bulence equations. Particularly, the additional terms in 
the Navier-Stokes equations due to the averaging pro- 
cedure are illustrated. These terms are often omitted 
in CFD practices. The k and c equations are presented 
in vector form for the convenience of illustrating the 
numerical method. 

The lower- upper (LU) scheme developed by Yoon 
and Jameson 2 is adopted in this work. This method has 
proven very efficient for large systems of equations. 3 For 
completeness, a brief account of the numerical method 
is presented in this paper. The newly developed solver 
then is applied to simulate compressible free shear lay- 
ers with five different convective Mach numbers from 
0.05 to 1.48. One of five test conditions is a replica 
of that reported by Goebel and Dutton. 4 The most 
important feature of the compressible free shear lay- 
ers one want to demonstrate in the calculations is the 
decrease of the shear layer thickness with increasing 
Mach number. 5 Favorable comparison were found be- 
tween the experimental data and the calculated results. 


THEORETICAL MODEL 

In deriving the compressible flow equations of fully 
turbulent flows, all the flow properties are Favre aver- 
aged (mass weighted averaged) except the density, p, 
and pressure, p. The definition of the Favre average is 

4> = 4>p/p • (1) 
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where is any flow properly. Thus flow variables are 
decomposed in the following fashion, 


= + 


( 2 ) 


On the other hand, conventional Reynolds’ average are 
used for the pressure and the density. According to 
the definition of Favre averaging, the following relations 
exist: 


r = 


— j 

p 


P<t > ” = o, 


(3) 


<y = o. 


In doing so, all the terms associated with the density 
fluctuation, e.g., p'u\ in the Reynolds’ averaged equa- 
tions were eliminated in the Favre ’s averaged equations. 
The resulting flow equations are much simpler com- 
pared to the equations derived by the Reynolds aver- 
aging procedures. 

Written in a strong conservative form, the turbu- 
lent, compressible, flow equations, can be expressed as 
follows: 


IT + |(E- E „) + |( F - F „)= H . 


(4) 


Here x and y are Cartesian coordinates, Q is the depen- 
dent variable, E and F are the convective flux vectors 
and E v and F v are the viscous flux vectors. The equa- 
tions are similar to the laminar equations. However, 
all variables in the equations are the averaged vari- 
ables, and the transport properties are the effective, 
i.e., laminar plus turbulent, properties. Here, we want 
to point out that the viscosity multiplied by the dilata- 
tional terms in the normal stresses is still the laminar 
viscosity as illustrated in the following relations: 


T XX 

Tyy 


0 , ,<9u 2 

+ --/i 

.dv 2 

= 2(v+p<)^--M 



(5) 


The vector H on the right hand side of the flow 
equations, Eq. 4, represents the additional terms intro- 
duced by the averaging procedure. The vector H can 
be expressed as 



V +2£t^(/>fc) + 2v^(pfc) + G / 


where k is the turbulent kinetic energy and is defined 
as 

k = ~rp(u”u” + (7) 

G is the generation of the turbulent kinetic energy and 
is defined in the following section. 

The equations of turbulent kinetic energy, Jb, and 
energy dissipation rate, €, are derived by the manipula- 
tion of the flow equations and the averaging procedure. 
The derived equations of k and c can not be solved di- 
rectly due to the closure problem. Modelling of certain 
terms in the equations is necessary to make the govern- 
ing equations well posed. The details of the modelling 
is beyond the scope of this work. Here, only the final 
form of the equations are presented. The k and c equa- 
tions in the Cartesian coordinate system can be cast 
into the vector form: . 


dQibe dEkc dF ke 

dt dx dy 


dE v kc , dJ? v kc 

dx dy 


+ S, 


( 8 ) 


where 



W+?.)£J’ 

(G-pe(l+aM?)\ 

V {C X G-Cipt)\ ) ’ 


^vke — 
F vkc ~ 


(9) 


where G is the generation term of the turbulent kinetic 
energy and can be expressed as: 


G = /i t 



~yi : <«» 


where the subscripts :, j, and k follows the convention 
of Cartesian tensor. 

The eddy viscosity p t is derived in terms of appro- 
priate length and velocity scales. For the k — e turbu- 
lence model, the length scale and the velocity scale of 
turbulent fluctuations are taken as 



« k l >\ 


( 11 ) 
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These relations allow the eddy viscosity, to be mod- 
elled as: 

( 12 ) 

The constants used in the k — c model are the standard 
Jones and Launder’s values : 6 C ^ = 0.09, C\ — 1.44, 
C 2 = 192, 4 X*. = 1.0, and cr t = 1.3. These constants 
were never altered during the course of this work. 

In order to accommodate the compressibility ef- 
fect, the dissipation term ~pc in the source term of the k 
equation is multiplied by a correction factor ( 1 -f aM ( 2 ). 
Here M t is the local turbulence Mach number defined 
as Mt = y/k/a where a is the local speed of sound. The 
constant oc in the term is taken as unity. This model 
is developed by Sarkar et al . 7 The physical meaning of 
the term is that for high turbulence Mach number (M t ) 
flows, the dissipation of the turbulence kinetic energy 
is enhanced by a factor of aM 2 . For free shear layers at 
high convective Mach number, the turbulence intensity 
is greatly reduced due to compressibility. 

In calculating the turbulent free shear layers, the 
inlet boundary conditions for mean velocities and tem- 
perature are specified based on the hyperbolic tangent 
profile with specified initial shear layer thickness. The 
hyperbolic tangent profile is an approximation of the 
self-similar solution for fully developed turbulent free 
shear layers. The inlet transverse velocities are set to 
be zero. The turbulent kinetic energy and dissipation 
are specified according to the local equilibrium assump- 
tion and a algebraic turbulence model : 8 

where / m = 0.1256 and the shear layer thickness, 6 , is 
based on the distance between the two transverse lo- 
cations where tt = tli — O.lAu and u = ^ -f O.lAti. 
The dissipation can be related to the local length scale 
which is specified based on the local shear layer thick- 
ness: 

c = C ( ^ (14) 

where C e = 1.23. Using Eqs. (13) and (14), and the 
equation for the eddy viscosity, i.e, Eq. (12), k and 
€ can be readily obtained for the upstream boundary 
conditions. 


Numerical Method 

The numerical solution of Eqs. (4) and ( 8 ) is 
performed in a general, body-fitted coordinate system, 
(£, 77 ). For the purpose of discussion, we will concen- 
trate on Eq. (4). However, the procedure is equally 


applicable to Eq. ( 8 ). Coordinate transformation of 
Eq. (4) results in: 

dQ . d ^ ^ \ . d 


+ 5 ?( E " Eu ) + ^( F-F '’) =I * (15) 


dr d£ 


where 


(16) 


dr] 

Q = /iQ 

E = /i(£ x E-KyF) 

F = h(r] x E + T] y F) 

= /i(^ r E v -f £yF„) 

Ft, = h(r] x E v + F v ) 

H = /iH 

in which h is the cell volume. 

The transformed equation, Eq. (15), is solved us- 
ing a time-marching, LU scheme. The LU scheme can 
be obtained by approximately factorizing the left-hand- 
side (LHS) of the equation. In time-marching form, the 
implicit upwind difference scheme of Eq. (15) can be 
written as 


[I - 1 - At(D^ A“ H- D+ B“ -D+ 

D~A+ -f D"B+)]AQ = AtRHS 


(17) 


In Eq. (17), At is the time-step. Backward-difference 
operators are denoted by D~ and D~, and forward- 
difference operators by and D + . The flux Jaco- 

bians, A + , B + , A“, and B“ are constructed such that 
the eigenvalues of matrices are nonnegative and 
those of matrices are nonpositive. The matrix D 
is the Jacobian matrix of the source term. 

The LHS matrix of Eq. (21) is usually too large for 
direct inversion. An approximate- factorization proce- 
dure is implemented which results in the following LU 
scheme : 


I + Al(o*A- + D+B- + ^-+|i-b 


AQ 


= AtR 

(18) 

where the grid spacing in the general coordinate, A£ 
and A 77 are usually taken to be one. R represents the 
residual of each LU time marching step. In calculating 
the residual, R, both the inviscid and viscous terms are 
discretized using the central-difference approach: 


R = Dz(E v - E) + D r] (F v - F) + H (19) 
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where D ^ and D n are the central difference operators. 

Equation (18) is the generic form for the LU 
scheme. Its derivation can be found in Ref. 3 and 
will not be repeated here. This LU scheme requires 
inversion of the matrix, 


is applied in this work to enhance the numerical stabil- 
ity. In the k equation, e has been replaced by k(e/k) 
where c/k is treated as a constant. A similar method is 
also applied to the source term of the e equation. The 
Jacobian matrix obtained is: 


[l + A< (A+-A- + B+ — B- — f))] (20) 


D = 


( 


-i(l + aM?) 
0 


° ) 

~C 


(27) 


for the L operator and 

[i + A«(a + -A- + B + -B--6/d)J (21) 


for the U operator. 

Up to this point, no definition has been made to 
the exact form of the split flux Jacobians. Yoon and 
Jameson 2 proposed that the split flux Jacobians are 
defined as 

A + =0.5(A + 7a I) 


A" = 0.5(A - 7aI) 
B+ = 0.5(B + 7 bI) 


( 22 ) 


B“ =0.5(B-76l) 


where 7^ and 73 are greater than the spectral radii of 
the associated flux Jacobians : 


7a > max (| A A|) 

7b > majc(|A 6 1) 


The purpose of constructing split flux-Jacobians by Eq. 
(22) is to make the matrices in Eqs. (20) and (21) diag- 
onal for efficient inversion. Apparently, the eigenvalues 
of the split flux-Jacobians are not the characteristics 
speeds of the flow. 

In solving the k and c equations, the aforemen- 
tioned numerical method, i.e., the LU scheme on the 
left hand side and central differencing on the right hand 
side, is used. The solution procedure of the whole equa- 
tion set is decoupled into flow solver and turbulence 
solver. Thus, the turbulence solver stands alone and 
can be easily turned on or off. This arrangement does 
not affect the overall numerical stability due to the fact 
that the feedback from k and e equations to the flow 
equations depends on the turbulent transport proper- 
ties only. Thus, it is more efficient and convenient to 
separate the solution procedure into two parts. 

The source terms of the k and e equations demand 
special treatment. In linearizing the source terms for 
the numerical method, the Jacobian matrix is obtained 
through the derivative of the source terms with respect 
to the dependent variables, i.e., ~pk and pe. Following 
the usual practice, the form of the source terms guar- 
antee a 2 x 2 full matrix for the Jacobian matrix. How- 
ever, special treatment in deriving the Jacobian matrix 


Note that off-diagonal terms are eliminated and the di- 
agonal terms are always negative. Thus, the implicit 
part of the source terms of k and e equations behaves 
like a sink which always stabilize the numerical scheme. 


Results and Discussions 

Bogdanoff 5 introduced the convective Mach num- 
ber as a parameter that collapses the growth rate date 
of plane shear layers. The convective Mach number is 
defined as 


Ui ~ Uc = U c ~ U 2 

Cl c 2 


(28) 


where U 1 and C/ 2 are the freestream velocities, and Ci 
and c 2 are the freestream sound speeds. U c is the con- 
vective velocity of large structures and is defined as 


U c = 


Uic 2 + U 2 C 1 

Cl + C 2 


(29) 


Table 1 shows the five test conditions of the simulated 
compressible free shear layers. The range of the convec- 
tive Mach numbers, M c , is from 0.05 to 1.48. The last 
two rows show the spreading rate and the ratio of the 
compressible spreading rate to incompressible spread- 
ing rate in which 6 is the shear layer thickness and the 
superscript ' represents the derivative of 6 with respect 
to the streamwise distance. As indicated in the last 
row of the table, the ratios of the spreading rate of 
the free shear layers decreases as the convective Mach 
number increases. The test conditions of Case 3 in the 
table are the same as in the experiments reported by 
Goebel and Dutton. 4 In the rest of the section, we will 
first show the direct comparison between the simulated 
results and the experimental data for Case 3. Then, 
detailed numerical solutions of Case 3 in terms of ve- 
locity, turbulent kinetic energy, turbulence dissipation 
rate, Reynolds stress, and eddy viscosity are presented 
in a coalesced fashion. Finally, the solutions of five 
cases are compared to each other in the figures of ki- 
netic energy, Reynolds stress, and ratio of spreading 
rates. 

The solutions of Case 3 are examined in detail. 
Figure 2 shows the development of the free shear layer. 
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Note that x and y axes are not on a 1:1 ratio for the 
convenience of illustration. The definition of the shear 
layer boundaries is the same that of the shear layer 
thickness. In Fig. 2, the boundaries of the shear layer 
corresponding to 10 — 90% are drawn. Circles are the 
experimental data of Dutton et al. The calculated re- 
sults agree well with the experimental data. After the 
developing region, the boundary of the shear layer is 
almost linear. Incidentally, Figs. 3a and 3b show the 
numerical convergence trends of the flow and k — c 
equations. In about 2500 iterations, the residuals drop 
about 12 orders of magnitude and reach the machine 
accuracy. 

Figure 4 shows the Mach number profile at various 
axial locations. The velocity gradient in the transverse 
direction decreases as the flow goes downstream. Figure 
5 is the coalesced version of Fig. 4. The nondimension- 
alized y coordinate defined as (y — y c )/6 is used, where 
y c is the transverse location at the center of the shear 
layer, and 6 is the shear layer thickness. Note that the 
upstream boundary condition of velocities is prescribed 
according to a hyperbolic tangent curve which is an ap- 
proximation of self-similar solution of free shear layers. 
According to Fig. 5, this self similarity of velocity pro- 
files never fail as the flow goes downstream. 

Figure 6 shows coalesced turbulence kinetic ener- 
gies at various locations. The turbulence kinetic en- 
ergy is nondimensionalized by AU 2 . This figure clearly 
shows that after the first three stations, i.e., about 150 
mm, the turbulence kinetic energy retains self similar- 
ity. Thus, the developing region for the turbulence is 
about 150 mm. Figure 7 shows the turbulence dissipa- 
tion profiles. As flow goes downstream, the peak values 
of turbulence dissipation at each stations decrease, and 
the turbulence dissipation never reach a fully developed 
condition. However, if the e value is nondimensional- 
ized by AU 3 /6 the coalesced profiles appear as shown 
in Fig. 7. A similar behavior is observed for the eddy 
viscosity profiles (Fig. 8). Figure 9 shows the nondi- 
mensionalized Reynolds stress profiles. Again, the tur- 
bulence Reynolds stress becomes the fully developed at 
about 150 mm downstream of the splitter plate. 

Figure 10 shows the comparison between the pre- 
dicted fully developed Reynolds stress and the experi- 
mental data reported by Dutton et al. 5 The predicted 
solution underestimated the peak value of the Reynolds 
stress profile by 6 ~ 8%; however, the overall trend 
of the predicted results is correct. Many factors con- 
tribute to the discrepancy between the predicted result 
and experimental data. Among them, the upstream 
boundary conditions are simplified in the solution pro- 
cedure, i.e., no effort was made to simulate two bound- 
ary layers merging at the tip of the splitter plate. This 
could offset the solution in the developing region and 


shift the fully developed solutions. 

Figure 11 shows the comparison of the turbulent 
kinetic energy between the five cases. Again, the turbu- 
lent kinetic energy is normalized by the square of the ve- 
locity difference of the two streams. It is clear that the 
normalized turbulence intensity decreases with increas- 
ing convective Mach number. A similar situation is 
observed in Fig. 12, the normalized Reynolds stress de- 
creases with increasing convective Mach number. Fig- 
ure 13 shows the distribution of the ratios of spread- 
ing rates for the five cases compared to experimental 
data. The ratio of the spreading rates decreases from 
about unity to 0.45 as the convective mach number in- 
creases from about 0. to 1.45. Both the experimental 
data and the simulated results show the spreading rate 
ratio reaches an asymptotic value after the convective 
mach number exceeds unity. Simular phenomenon can 
be seen in Figs. 11 and 12. Both figures show that 
the normalized turbulence kinetic energy and Reynolds 
stress reach asymptotic values as the convective Mach 
number increases. 


CONCLUDING REMARKS 

In this paper, we report the incorporation of a com- 
pressible k — c model with a two-dimensional Navier 
Stokes solver to study compressible free shear layers. 
Sarkar’s modelling is adopted to simulate the compress- 
ibility effects of the k and c equations. The model en- 
hances the turbulence dissipation rate of flows at high 
speeds. In deriving the governing equations, the Favre 
averaging procedure for the fully trubulent flow equa- 
tions is elaborated. The equation sets are presented in 
the vector form for the convenience of the discussion of 
the numerical method. Yoon and Jameson’s LU scheme 
is used to solve the equation sets. Details of boundary 
conditions and the treatment of the source terms of the 
k and c equations are discussed. Then the program is 
applied to simulate compressible free shear layers with 
five different convective Mach numbers. The decrease 
of the spreading rate of the shear layers with increasing 
Mach number is observed in the calculated results. Re- 
sults also show favorable comparison with Goebel and 
Dutton’s experimental data. 
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u 2 


Fig. 1: Schematic of free shear layer. 


Table T. Test conditions of the calculated shear layer. 


CASE 

1 

2 

3 

4 

5 

Mj,M 2 

1.2, 1.1 

2.0, 1.4 

1.91, 1.37 

2.5, 1.1 

6.1, 3.15 

M c 

0.05 

0.31 

0.45 

0.70 

1.48 

u 15 u 2 

[m] 

419, 384 

676, 465 

702, 404 

865, 380 

3461, 1786 

T * T * 

300, 300 

275, 275 

334, 215 

295, 295 

800, 800 

Pi* P? 

(Kg/m*3] 

0.64,0.6^ 

0.7,0. 7 

0.57,0.89 

0.64,0.64 

0.24,0.24 

P [atm] 

0.55 

0.55 

0.55 

0.55 

0.55 

8’ 

0.007 

0.021 

0.027 

0.035 

0.024 

876’ L 

1.01 

0.7 

0.54 

0.47 

0.40 
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1.1 I 


Fig. 4 : The velocity profiles of the free shear layer at dif- 
ferent axial locations. 


Fig. 2: The boundaries of the free shear layer (10 90%). 
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Fig. £: Coalesced velocity 
tions. 
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o 348 . 33 
a 398.33 
v 448 . 33 
► 498.33 

profiles of different axial loca- 


Fig. 3: Convergence trends of the flow and turbulence 
equations. 
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Fig. 6: Coalesced turbulent kinetic energy profiles of dif- 
ferent axial locations. 



Fig. 7: Coalesced turbulence dissipation profiles of differ- 
ent axial locations. 


Fig. 9: Coalesced Reynolds stress profiles of different axial 
locations. 



Fig. 10: Comparison of Reynolds stresses between experi- 
mental data and predicted results. 



Fig. 8: Coalesced eddy viscosity profiles of different axial 
locations. 
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ig. 11: Comparison of turbulent kinetic energy with five 
different convective Mach numbers. 



A Papamoschou and Roshko [9] 

Q Langley’s data [7] 

-Q— With Compressibility Model 

Without Compressibility Model 

Fig. 13: Comparison of ratios of spreading rate between ex- 
perimental data and predicted results. 


Fig. 12: Comparison of Reynolds stresses with five different 
convective Mach numbers. 
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